AtomSet.C 16.7 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19
// AtomSet.C
//
////////////////////////////////////////////////////////////////////////////////

#include "AtomSet.h"
20
#include "Species.h"
Francois Gygi committed
21
#include "NameOf.h"
22
#include "sampling.h"
Francois Gygi committed
23 24 25 26 27
#include <iostream>
#include <algorithm>
#include <string>
using namespace std;

28 29 30 31 32 33 34 35 36 37 38
////////////////////////////////////////////////////////////////////////////////
AtomSet::~AtomSet(void)
{
  for ( int is = 0; is < species_list.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      delete atom_list[is][ia];
    }
    delete species_list[is];
  }
Francois Gygi committed
39
}
Francois Gygi committed
40 41 42
////////////////////////////////////////////////////////////////////////////////
bool AtomSet::addSpecies(Species* sp, string name)
{
Francois Gygi committed
43 44 45 46 47
  const double rcps = 1.5;
  sp->initialize(rcps);

  Species *s = findSpecies(name);
  if ( s != 0 )
Francois Gygi committed
48
  {
Francois Gygi committed
49
    // species is already defined: substitute with new definition
Francois Gygi committed
50
    if ( ctxt_.onpe0() )
Francois Gygi committed
51
    {
Francois Gygi committed
52 53
      cout << " AtomSet::addSpecies: species " << name
           << " is already defined" << endl;
Francois Gygi committed
54 55 56 57 58 59 60 61 62
      cout << " AtomSet::addSpecies: redefining species" << endl;
    }
    // Check if s and sp are compatible: zval must be equal
    if ( s->zval() != sp->zval() )
    {
      cout << " AtomSet::addSpecies: species do not have the same"
           << " number of valence electrons. Cannot redefine species" << endl;
      return false;
    }
63

Francois Gygi committed
64 65 66 67 68 69 70 71 72 73 74 75 76
    int is = isp_[s->name()];
    species_list[is] = sp;
    delete s;
  }
  else
  {
    // create new entry in species list
    species_list.push_back(sp);
    spname.push_back(name);
    isp_[name] = species_list.size()-1;
    na_[name] = 0;
    atom_list.resize(atom_list.size()+1);
  }
77

78 79
  if ( ctxt_.onpe0() )
  {
80
    cout << endl << " species " << sp->name() << ":" << endl;
81 82
    sp->info(cout);
  }
Francois Gygi committed
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
  return true;
}

////////////////////////////////////////////////////////////////////////////////
bool AtomSet::addAtom(Atom *a)
{

  // check atom_list for name
  if ( findAtom(a->name()) )
  {
    // this name is already in the atom list, reject atom definition
    if ( ctxt_.onpe0() )
      cout << " AtomSet:addAtom: atom " << a->name()
           << " is already defined" << endl;
    return false;
  }

  // check if species is defined
  string spname = a->species();
  Species *s = findSpecies(spname);
  if ( !s )
  {
    // species not found, cannot define atom
    if ( ctxt_.onpe0() )
      cout << " AtomSet:addAtom: species " << spname
           << " is undefined" << endl;
    return false;
  }
111

Francois Gygi committed
112 113 114 115
  // add an atom to the atom_list
  int is = isp(spname);
  assert ( is >= 0 );
  atom_list[is].push_back(a);
116 117
  ia_[a->name()] = atom_list[is].size()-1;
  is_[a->name()] = is;
Francois Gygi committed
118 119 120 121

  // update count of atoms of species spname
  // increment count
  na_[spname]++;
122

Francois Gygi committed
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
  // update total number of electrons
  nel_ = 0;
  for ( int is = 0; is < species_list.size(); is++ )
  {
    nel_ += species_list[is]->zval() * atom_list[is].size();
  }

  return true;
}

////////////////////////////////////////////////////////////////////////////////
bool AtomSet::delAtom(string name)
{
  vector<vector<Atom*> >::iterator psa = atom_list.begin();
  vector<Species*>::iterator ps = species_list.begin();
  for ( int is = 0; is < species_list.size(); is++ )
  {
    vector<Atom*>::iterator pa = atom_list[is].begin();
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      if ( atom_list[is][ia]->name() == name )
      {
        string spname = atom_list[is][ia]->species();
        na_[spname]--;
        delete atom_list[is][ia];
148

149 150 151 152 153
        // remove map entries ia_[name] and is_[name]
        map<string,int>::iterator i = ia_.find(name);
        ia_.erase(i);
        i = is_.find(name);
        is_.erase(i);
154

Francois Gygi committed
155 156
        atom_list[is].erase(pa);
        nel_ -= species_list[is]->zval();
157

Francois Gygi committed
158 159 160 161 162 163 164
        return true;
      }
      pa++;
    }
    psa++;
    ps++;
  }
165

Francois Gygi committed
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
  // this name was not found in the atom list
  if ( ctxt_.onpe0() )
    cout << " AtomSet:delAtom: no such atom: " << name << endl;
  return false;
}

////////////////////////////////////////////////////////////////////////////////
Atom *AtomSet::findAtom(string name) const
{
  for ( int is = 0; is < species_list.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      if ( atom_list[is][ia]->name() == name )
        return atom_list[is][ia];
    }
  }
  return 0;
}

////////////////////////////////////////////////////////////////////////////////
Species *AtomSet::findSpecies(string name) const
{
  int is = isp(name);
  if ( is >= 0 )
    return species_list[is];
  else
    return 0;
}

////////////////////////////////////////////////////////////////////////////////
197
int AtomSet::isp(const string& name) const
Francois Gygi committed
198 199 200 201 202 203 204 205
{
  map<string,int>::const_iterator i = isp_.find(name);
  if ( i != isp_.end() )
    return (*i).second;
  else
    return -1;
}

206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
////////////////////////////////////////////////////////////////////////////////
int AtomSet::is(const string& atom_name) const
{
  map<string,int>::const_iterator i = is_.find(atom_name);
  if ( i != is_.end() )
    return (*i).second;
  else
    return -1;
}

////////////////////////////////////////////////////////////////////////////////
int AtomSet::ia(const string& atom_name) const
{
  map<string,int>::const_iterator i = ia_.find(atom_name);
  if ( i != ia_.end() )
    return (*i).second;
  else
    return -1;
}

Francois Gygi committed
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
////////////////////////////////////////////////////////////////////////////////
void AtomSet::listAtoms(void) const
{
  for ( int is = 0; is < species_list.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      if ( ctxt_.onpe0() )
        cout << *atom_list[is][ia];
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void AtomSet::listSpecies(void) const
{
  for ( int is = 0; is < species_list.size(); is++ )
  {
    if ( ctxt_.onpe0() )
    {
      cout << endl << " species " << spname[is] << ":" << endl;
      species_list[is]->info(cout);
      cout << endl;
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
254
int AtomSet::na(const string& spname) const
Francois Gygi committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
{
  map<string,int>::const_iterator i = na_.find(spname);
  if ( i != na_.end() )
    return (*i).second;
  else
    return 0;
}

////////////////////////////////////////////////////////////////////////////////
int AtomSet::na(int is) const
{
  assert( is >= 0 && is < atom_list.size() );
  return atom_list[is].size();
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
271
int AtomSet::size(void) const
Francois Gygi committed
272 273 274 275 276 277 278 279 280 281 282
{
  int n = 0;
  for ( int is = 0; is < atom_list.size(); is++ )
    n += atom_list[is].size();
  return n;
}

////////////////////////////////////////////////////////////////////////////////
bool AtomSet::reset(void)
{
  // delete all atoms and species
283

Francois Gygi committed
284 285 286 287 288 289 290
  for ( int is = 0; is < species_list.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      delete atom_list[is][ia];
    }
    atom_list[is].resize(0);
291

Francois Gygi committed
292 293 294 295
    delete species_list[is];
  }
  atom_list.resize(0);
  species_list.resize(0);
Francois Gygi committed
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
  spname.resize(0);
  for ( map<string,int>::iterator i = na_.begin();
        i != na_.end();
        i++ )
        na_.erase(i);

  for ( map<string,int>::iterator i = isp_.begin();
        i != isp_.end();
        i++ )
        isp_.erase(i);

  for ( map<string,int>::iterator i = is_.begin();
        i != is_.end();
        i++ )
        is_.erase(i);
311

Francois Gygi committed
312 313 314 315
  for ( map<string,int>::iterator i = ia_.begin();
        i != ia_.end();
        i++ )
        ia_.erase(i);
Francois Gygi committed
316
  nel_ = 0;
317

Francois Gygi committed
318 319 320 321 322 323
  return true;
}

////////////////////////////////////////////////////////////////////////////////
void AtomSet::get_positions(vector<vector<double> >& tau) const
{
324 325
  if (tau.size() != atom_list.size())
    tau.resize(atom_list.size());
Francois Gygi committed
326 327
  for ( int is = 0; is < atom_list.size(); is++ )
  {
328 329
    if (tau[is].size() != 3*atom_list[is].size())
      tau[is].resize(3*atom_list[is].size());
Francois Gygi committed
330 331 332 333 334 335 336 337 338 339 340
    int i = 0;
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector t = atom_list[is][ia]->position();
      tau[is][i++] = t.x;
      tau[is][i++] = t.y;
      tau[is][i++] = t.z;
    }
  }
}

Francois Gygi committed
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358
////////////////////////////////////////////////////////////////////////////////
void AtomSet::sync_positions(vector<vector<double> >& tau)
{
  for ( int is = 0; is < tau.size(); is++ )
  {
    int m = tau[is].size();
    double* p = &tau[is][0];
    if ( ctxt_.onpe0() )
    {
      ctxt_.dbcast_send(m,1,p,m);
    }
    else
    {
      ctxt_.dbcast_recv(m,1,p,m,0,0);
    }
  }
}

Francois Gygi committed
359
////////////////////////////////////////////////////////////////////////////////
360
void AtomSet::set_positions(vector<vector<double> >& tau)
Francois Gygi committed
361
{
362
  sync_positions(tau);
Francois Gygi committed
363 364 365 366
  assert(tau.size() == atom_list.size());
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    assert(tau[is].size() == 3*atom_list[is].size());
Francois Gygi committed
367
    for ( int ia = 0, i = 0; ia < atom_list[is].size(); ia++ )
Francois Gygi committed
368 369 370 371 372 373 374 375 376 377 378
    {
      atom_list[is][ia]->set_position(
        D3vector(tau[is][i],tau[is][i+1],tau[is][i+2]));
      i += 3;
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void AtomSet::get_velocities(vector<vector<double> >& vel) const
{
379 380
  if (vel.size() != atom_list.size())
    vel.resize(atom_list.size());
Francois Gygi committed
381 382
  for ( int is = 0; is < atom_list.size(); is++ )
  {
383 384
    if (vel[is].size() != 3*atom_list[is].size())
      vel[is].resize(3*atom_list[is].size());
Francois Gygi committed
385 386 387 388 389 390 391 392 393 394 395
    int i = 0;
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector t = atom_list[is][ia]->velocity();
      vel[is][i++] = t.x;
      vel[is][i++] = t.y;
      vel[is][i++] = t.z;
    }
  }
}

Francois Gygi committed
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
////////////////////////////////////////////////////////////////////////////////
void AtomSet::sync_velocities(vector<vector<double> >& vel)
{
  for ( int is = 0; is < vel.size(); is++ )
  {
    int m = vel[is].size();
    double* p = &vel[is][0];
    if ( ctxt_.onpe0() )
    {
      ctxt_.dbcast_send(m,1,p,m);
    }
    else
    {
      ctxt_.dbcast_recv(m,1,p,m,0,0);
    }
  }
}

Francois Gygi committed
414
////////////////////////////////////////////////////////////////////////////////
415
void AtomSet::set_velocities(vector<vector<double> >& vel)
Francois Gygi committed
416
{
417
  sync_velocities(vel);
Francois Gygi committed
418 419 420 421
  assert(vel.size() == atom_list.size());
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    assert(vel[is].size() == 3*atom_list[is].size());
Francois Gygi committed
422
    for ( int ia = 0, i = 0; ia < atom_list[is].size(); ia++ )
Francois Gygi committed
423 424 425 426 427 428 429 430
    {
      atom_list[is][ia]->set_velocity(
        D3vector(vel[is][i],vel[is][i+1],vel[is][i+2]));
      i += 3;
    }
  }
}

431 432 433 434 435 436 437 438 439 440
////////////////////////////////////////////////////////////////////////////////
void AtomSet::reset_velocities(void)
{
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
      atom_list[is][ia]->set_velocity(D3vector(0.0, 0.0, 0.0));
  }
}

441 442 443 444 445 446 447
////////////////////////////////////////////////////////////////////////////////
void AtomSet::rescale_velocities(double fac)
{
  vector<vector<double> > v;
  get_velocities(v);
  for ( int is = 0; is < v.size(); is++ )
  {
448 449
    for ( int i = 0; i < v[is].size(); i++ )
      v[is][i] *= fac;
450 451 452 453
  }
  set_velocities(v);
}

454 455 456
////////////////////////////////////////////////////////////////////////////////
void AtomSet::randomize_positions(double amplitude)
{
Francois Gygi committed
457
  // add random displacements to positions using
458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
  // random numbers from a normal distribution scaled
  // by the amplitude parameter
  vector<vector<double> > r;
  get_positions(r);
  for ( int is = 0; is < r.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      // draw pairs of unit variance gaussian random variables
      double xi0, xi1, xi2, xi3; // xi3 not used
      normal_dev(&xi0,&xi1);
      normal_dev(&xi2,&xi3);
      r[is][3*ia+0] += amplitude * xi0;
      r[is][3*ia+1] += amplitude * xi1;
      r[is][3*ia+2] += amplitude * xi2;
    }
  }
  set_positions(r);
}

478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501
////////////////////////////////////////////////////////////////////////////////
void AtomSet::randomize_velocities(double temp)
{
  // initialize velocities with random numbers from a Maxwell-Boltzmann
  // distribution at temperature temp
  const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
  vector<vector<double> > v;
  get_velocities(v);
  for ( int is = 0; is < v.size(); is++ )
  {
    const double m = species_list[is]->mass() * 1822.89;
    const double width = sqrt( boltz * temp / m );
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      // draw pairs of unit variance gaussian random variables
      double xi0, xi1, xi2, xi3; // xi3 not used
      normal_dev(&xi0,&xi1);
      normal_dev(&xi2,&xi3);
      v[is][3*ia+0] = width * xi0;
      v[is][3*ia+1] = width * xi1;
      v[is][3*ia+2] = width * xi2;
    }
  }
  set_velocities(v);
502
  reset_vcm();
503
}
504

505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::vcm(void) const
{
  D3vector mvsum;
  double msum = 0.0;
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    double mass = species_list[is]->mass();
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector v = atom_list[is][ia]->velocity();
      mvsum += mass * v;
      msum += mass;
    }
  }
  if ( msum == 0.0 ) return D3vector(0,0,0);
  return mvsum / msum;
}

////////////////////////////////////////////////////////////////////////////////
void AtomSet::reset_vcm(void)
{
  D3vector vc = vcm();
  vector<vector<double> > v;
  get_velocities(v);
  // subtract center of mass velocity
  for ( int is = 0; is < v.size(); is++ )
  {
    int i = 0;
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      v[is][i++] -= vc.x;
      v[is][i++] -= vc.y;
      v[is][i++] -= vc.z;
    }
  }
  set_velocities(v);
}

Francois Gygi committed
544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562
////////////////////////////////////////////////////////////////////////////////
void AtomSet::fold_in_ws(void)
{
  vector<vector<double> > p;
  get_positions(p);
  for ( int is = 0; is < p.size(); is++ )
  {
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector pos(&p[is][3*ia]);
      cell_.fold_in_ws(pos);
      p[is][3*ia+0] = pos.x;
      p[is][3*ia+1] = pos.y;
      p[is][3*ia+2] = pos.z;
    }
  }
  set_positions(p);
}

Francois Gygi committed
563 564 565 566 567 568 569 570 571 572
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::dipole(void) const
{
  D3vector sum;
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    double charge = species_list[is]->zval();
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector p = atom_list[is][ia]->position();
573
      sum += charge * p;
Francois Gygi committed
574 575 576 577 578
    }
  }
  return sum;
}

579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596
////////////////////////////////////////////////////////////////////////////////
D3tensor AtomSet::quadrupole(void) const
{
  D3tensor sum;
  for ( int is = 0; is < atom_list.size(); is++ )
  {
    double charge = species_list[is]->zval();
    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
    {
      D3vector p = atom_list[is][ia]->position();
      for ( int idir = 0; idir < 3; idir++ )
        for ( int jdir = 0; jdir < 3; jdir++ )
          sum[idir*3+jdir] += charge * p[idir] * p[jdir];
    }
  }
  return sum;
}

597
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
598
void AtomSet::set_cell(const D3vector& a, const D3vector& b, const D3vector& c)
599
{
Francois Gygi committed
600
  cell_.set(a,b,c);
601
  sync_cell();
Francois Gygi committed
602 603 604 605 606 607 608
}

////////////////////////////////////////////////////////////////////////////////
void AtomSet::sync_cell(void)
{
  sync_cell(cell_);
}
609

Francois Gygi committed
610 611 612
////////////////////////////////////////////////////////////////////////////////
void AtomSet::sync_cell(UnitCell& cell)
{
613 614 615 616
  if ( ctxt_.onpe0() )
  {
    double sbuf[9];
    for ( int i = 0; i < 9; i++ )
Francois Gygi committed
617
      sbuf[i] = cell.amat(i);
618 619 620 621 622 623 624 625 626
    ctxt_.dbcast_send(9,1,sbuf,9);
  }
  else
  {
    double rbuf[9];
    ctxt_.dbcast_recv(9,1,rbuf,9,0,0);
    D3vector a0(rbuf[0],rbuf[1],rbuf[2]);
    D3vector a1(rbuf[3],rbuf[4],rbuf[5]);
    D3vector a2(rbuf[6],rbuf[7],rbuf[8]);
Francois Gygi committed
627
    cell.set(a0,a1,a2);
628
  }
629
}
Francois Gygi committed
630

Francois Gygi committed
631
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
632
ostream& operator << ( ostream &os, const AtomSet &as )
Francois Gygi committed
633 634 635 636
{
  if ( as.context().onpe0() )
  {
    os << "<atomset>\n";
Francois Gygi committed
637
    os << as.cell();
Francois Gygi committed
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653
    for ( int is = 0; is < as.species_list.size(); is++ )
    {
      os << *as.species_list[is];
    }
    for ( int is = 0; is < as.species_list.size(); is++ )
    {
      for ( int ia = 0; ia < as.atom_list[is].size(); ia++ )
      {
        os << *as.atom_list[is][ia];
      }
    }
  }
  os << "</atomset>\n";

  return os;
}