Commit c6074490 by Francois Gygi

rel1_51_0


git-svn-id: http://qboxcode.org/svn/qb/trunk@764 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 685ad1c1
......@@ -15,7 +15,7 @@
// AtomSet.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomSet.C,v 1.27 2009-10-06 06:30:27 fgygi Exp $
// $Id: AtomSet.C,v 1.28 2010-02-20 23:13:02 fgygi Exp $
#include "AtomSet.h"
#include "Species.h"
......@@ -40,23 +40,40 @@ AtomSet::~AtomSet(void)
////////////////////////////////////////////////////////////////////////////////
bool AtomSet::addSpecies(Species* sp, string name)
{
if ( findSpecies(name) )
const double rcps = 1.5;
sp->initialize(rcps);
Species *s = findSpecies(name);
if ( s != 0 )
{
// species is already defined: substitute with new definition
if ( ctxt_.onpe0() )
{
cout << " AtomSet::addSpecies: species " << name
<< " is already defined" << endl;
return false;
}
const double rcps = 1.5;
sp->initialize(rcps);
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;
}
// 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);
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);
}
if ( ctxt_.onpe0() )
{
......@@ -276,7 +293,26 @@ bool AtomSet::reset(void)
}
atom_list.resize(0);
species_list.resize(0);
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);
for ( map<string,int>::iterator i = ia_.begin();
i != ia_.end();
i++ )
ia_.erase(i);
nel_ = 0;
return true;
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// AtomicExtForce.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomicExtForce.C,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#include "AtomicExtForce.h"
#include "AtomSet.h"
#include "Atom.h"
#include "Species.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void AtomicExtForce::setup(const AtomSet& atoms)
{
// find position in tau array corresponding to atom name1
is1_ = atoms.is(name1_);
ia1_ = atoms.ia(name1_);
assert(is1_>=0);
assert(ia1_>=0);
}
////////////////////////////////////////////////////////////////////////////////
double AtomicExtForce::energy(const vector<vector<double> > &r,
vector<vector<double> > &f)
{
f[is1_][3*ia1_+0] += force_.x;
f[is1_][3*ia1_+1] += force_.y;
f[is1_][3*ia1_+2] += force_.z;
return -force_.x*r[is1_][3*ia1_+0] +
force_.y*r[is1_][3*ia1_+1] +
force_.z*r[is1_][3*ia1_+2];
}
////////////////////////////////////////////////////////////////////////////////
ostream& AtomicExtForce::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <extforce name=\"" << name();
os << "\" type=\"" << type();
os << "\" atoms=\"" << name1_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << force_.x << " "
<< setprecision(6) << force_.y << " "
<< setprecision(6) << force_.z << "\"/>";
return os;
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// AtomicExtForce.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomicExtForce.h,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#ifndef ATOMICEXTFORCE_H
#define ATOMICEXTFORCE_H
#include "ExtForce.h"
#include "D3vector.h"
#include <cassert>
class AtomSet;
class AtomicExtForce: public ExtForce
{
std::string name1_; // name of atom
int ia1_, is1_;
D3vector force_;
public:
AtomicExtForce(std::string name, std::string name1, D3vector force):
name1_(name1), force_(force)
{
name_ = name;
names_.resize(1);
names_[0] = name1_;
}
~AtomicExtForce(void) {}
std::string type(void) const { return "atomic"; }
void set_value(std::vector<double> f)
{
force_.x = f[0];
force_.y = f[1];
force_.z = f[2];
}
void setup(const AtomSet& atoms);
double energy(const std::vector<std::vector<double> > &r,
std::vector<std::vector<double> > &f);
D3vector sum_contrib(void) { return force_; }
double magnitude(void) { return length(force_); }
std::ostream& print( std::ostream& os );
};
#endif
......@@ -15,7 +15,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.54 2009-09-08 18:45:59 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.55 2010-02-20 23:13:02 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -389,8 +389,10 @@ void BOSampleStepper::step(int niter)
<< " <exc> " << setw(15) << ef_.exc() << " </exc>\n"
<< " <esr> " << setw(15) << ef_.esr() << " </esr>\n"
<< " <eself> " << setw(15) << ef_.eself() << " </eself>\n"
<< " <ets> " << setw(15) << ef_.ets() << " </ets>\n"
<< " <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n";
<< " <ets> " << setw(15) << ef_.ets() << " </ets>\n";
if ( s_.extforces.size() > 0 )
cout << " <eexf> " << setw(15) << ef_.eexf() << " </eexf>\n";
cout << " <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n";
if ( compute_stress )
{
const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
......
......@@ -15,7 +15,7 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.22 2009-05-15 04:37:02 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.23 2010-02-20 23:13:02 fgygi Exp $
#include "CPSampleStepper.h"
#include "SlaterDet.h"
......@@ -166,8 +166,10 @@ void CPSampleStepper::step(int niter)
<< " <ecoul> " << setw(15) << ef_.ecoul() << " </ecoul>\n"
<< " <exc> " << setw(15) << ef_.exc() << " </exc>\n"
<< " <esr> " << setw(15) << ef_.esr() << " </esr>\n"
<< " <eself> " << setw(15) << ef_.eself() << " </eself>\n"
<< " <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n"
<< " <eself> " << setw(15) << ef_.eself() << " </eself>\n";
if ( s_.extforces.size() > 0 )
cout << " <eexf> " << setw(15) << ef_.eexf() << " </eexf>\n";
cout << " <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n"
<< flush;
if ( compute_stress )
{
......
......@@ -15,7 +15,7 @@
// ConstraintSet.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.C,v 1.10 2009-10-06 06:23:28 fgygi Exp $
// $Id: ConstraintSet.C,v 1.11 2010-02-20 23:13:02 fgygi Exp $
#include "ConstraintSet.h"
#include "PositionConstraint.h"
......@@ -643,3 +643,12 @@ void ConstraintSet::setup(AtomSet& atoms)
constraint_list[i]->setup(atoms);
}
}
////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::reset(void)
{
for ( int i = 0; i < constraint_list.size(); i++ )
delete constraint_list[i];
ndofs_ = 0;
constraint_list.resize(0);
}
......@@ -15,7 +15,7 @@
// ConstraintSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.h,v 1.9 2010-01-09 19:35:40 fgygi Exp $
// $Id: ConstraintSet.h,v 1.10 2010-02-20 23:13:02 fgygi Exp $
#ifndef CONSTRAINTSET_H
#define CONSTRAINTSET_H
......@@ -56,5 +56,6 @@ class ConstraintSet
const std::vector<std::vector<double> > &f);
void update_constraints(double dt);
void setup(AtomSet& atoms);
void reset(void);
};
#endif
......@@ -15,7 +15,7 @@
// EnergyFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.C,v 1.35 2009-11-05 07:03:39 fgygi Exp $
// $Id: EnergyFunctional.C,v 1.36 2010-02-20 23:13:02 fgygi Exp $
#include "EnergyFunctional.h"
#include "Sample.h"
......@@ -118,6 +118,7 @@ EnergyFunctional::EnergyFunctional(const Sample& s, const ChargeDensity& cd)
tau0.resize(nsp_);
fion_esr.resize(nsp_);
fext.resize(nsp_);
ftmp.resize(3*namax_);
eself_ = 0.0;
......@@ -129,6 +130,7 @@ EnergyFunctional::EnergyFunctional(const Sample& s, const ChargeDensity& cd)
const int na = atoms.na(is);
tau0[is].resize(3*na);
fion_esr[is].resize(3*na);
fext[is].resize(3*na);
eself_ += na * s->eself();
na_[is] = na;
......@@ -161,7 +163,6 @@ EnergyFunctional::EnergyFunctional(const Sample& s, const ChargeDensity& cd)
cell_moved();
atoms_moved();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -314,8 +315,8 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const bool use_confinement = s_.ctrl.ecuts > 0.0;
for ( int is = 0; is < fion.size(); is++ )
for ( int ia = 0; ia < fion[is].size(); ia++ )
fion[is][ia] = 0.0;
for ( int i = 0; i < fion[is].size(); i++ )
fion[is][i] = 0.0;
if ( compute_hpsi )
{
......@@ -603,7 +604,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
ets_ = - wf_entropy * s_.ctrl.fermi_temp * boltz;
}
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_;
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
if ( compute_hpsi )
{
......@@ -726,6 +727,15 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
{
fion[is][i] += fion_esr[is][i] + fac * ftmp[i];
}
// add external forces
if ( s_.extforces.size() > 0 )
{
assert(fion.size()==fext.size());
assert(fion[is].size()==fext[is].size());
for ( int i = 0; i < fion[is].size(); i++ )
fion[is][i] += fext[is][i];
}
}
}
......@@ -839,7 +849,6 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
<< " <sigma_eks_xz> " << setw(12) << sigma[5] << " </sigma_eks_xz>\n"
<< " </stress_tensor>" << endl;
}
return etotal_;
}
......@@ -967,6 +976,10 @@ void EnergyFunctional::atoms_moved(void)
}
}
sigma_esr *= - omega_inv;
// get external forces in fext
eexf_ = s_.extforces.energy(tau0,fext);
}
////////////////////////////////////////////////////////////////////////////////
......@@ -1020,6 +1033,7 @@ void EnergyFunctional::print(ostream& os) const
<< " <esr> " << setw(15) << esr() << " </esr>\n"
<< " <eself> " << setw(15) << eself() << " </eself>\n"
<< " <ets> " << setw(15) << ets() << " </ets>\n"
<< " <eexf> " << setw(15) << eexf() << " </eexf>\n"
<< " <etotal> " << setw(15) << etotal() << " </etotal>\n"
<< flush;
}
......
......@@ -15,7 +15,7 @@
// EnergyFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.h,v 1.19 2008-09-08 15:56:18 fgygi Exp $
// $Id: EnergyFunctional.h,v 1.20 2010-02-20 23:13:02 fgygi Exp $
#ifndef ENERGYFUNCTIONAL_H
#define ENERGYFUNCTIONAL_H
......@@ -61,12 +61,13 @@ class EnergyFunctional
std::vector<double> ftmp;
std::vector<std::vector<double> > tau0, fion_esr;
std::vector<std::vector<double> > fext;
std::vector<double> zv_, rcps_;
std::vector<int> na_;
int namax_;
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
ecoul_, exc_, esr_, eself_, ets_, etotal_;
ecoul_, exc_, esr_, eself_, ets_, eexf_, etotal_;
std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
......@@ -90,6 +91,7 @@ class EnergyFunctional
double esr(void) const { return esr_; }
double eself(void) const { return eself_; }
double ets(void) const { return ets_; }
double eexf(void) const { return eexf_; }
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ExtForce.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ExtForce.C,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#include "ExtForce.h"
#include <iostream>
using namespace std;
ostream& operator << ( ostream &os, ExtForce &x )
{
return x.print(os);
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ExtForce.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ExtForce.h,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#ifndef EXTFORCE_H
#define EXTFORCE_H
#include <string>
#include <vector>
#include <cassert>
#include "D3vector.h"
class AtomSet;
class ExtForce
{
protected:
std::string name_; // extforce name
std::vector<std::string> names_; // names of atoms involved in the extforce
public:
virtual ~ExtForce(void){}
virtual std::string type(void) const = 0;
virtual void set_value(std::vector<double> f) = 0;
virtual double energy(const std::vector<std::vector<double> > &r,
std::vector<std::vector<double> > &f) = 0;
virtual void setup(const AtomSet& atoms) = 0;
virtual std::ostream& print(std::ostream &os) = 0;
virtual D3vector sum_contrib(void) = 0;
std::string name(void) const { return name_; }
virtual double magnitude(void) = 0;
std::string names(int i) const
{
assert( i >= 0 && i < names_.size() );
return names_[i];
}
};
std::ostream& operator << ( std::ostream &os, ExtForce &xf );
#endif
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ExtForceCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ExtForceCmd.h,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#ifndef EXTFORCECMD_H
#define EXTFORCECMD_H
#include "UserInterface.h"
#include "Sample.h"
class ExtForceCmd : public Cmd
{
public:
Sample *s;
ExtForceCmd(Sample *sample) : s(sample) {};
char *name(void) const { return "extforce"; }
char *help_msg(void) const
{
return
"\n extforce\n\n"
" syntax:\n\n"
" extforce define atomic name atom fx fy fz\n"
" extforce define pair name atom1 atom2 force\n"
" extforce define global name fx fy fz\n"
" extforce set name value fx fy fz\n"
" extforce set name value f\n"
" extforce delete name\n"
" extforce list\n"
" External forces are added to ionic forces at each MD step.\n\n";
}
int action(int argc, char **argv)
{
const bool onpe0 = s->ctxt_.onpe0();
if ( argc < 2 )
{
if ( onpe0 )
cout << help_msg();
return 1;
}
string subcmd(argv[1]);
if ( subcmd == "define" )
{
return s->extforces.define_extforce(s->atoms,argc,argv);
}
else if ( subcmd == "set" )
{
return s->extforces.set_extforce(argc,argv);
}
else if ( subcmd == "delete" )
{
return s->extforces.delete_extforce(argc,argv);
}
else if ( subcmd == "list" )
{
if ( onpe0 )
s->extforces.list_extforces(cout);
}
else
{
if ( onpe0 )
cout << help_msg();
}
return 0;
}
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ExtForceSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ExtForceSet.h,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#ifndef EXTFORCESET_H
#define EXTFORCESET_H
#include <vector>
#include <string>
class Atom;
class AtomSet;
class ExtForce;
class Context;
class ExtForceSet
{
private:
const Context& ctxt_;
std::vector<ExtForce *> extforce_list;
public:
ExtForceSet(const Context& ctxt) : ctxt_(ctxt) {}
~ExtForceSet();
bool define_extforce(AtomSet &atoms, int argc, char **argv);
bool set_extforce(int argc, char **argv);
bool delete_extforce(int argc, char **argv);
void list_extforces(std::ostream &os);
int size(void) const { return extforce_list.size(); }
double energy(const std::vector<std::vector<double> > &r,
std::vector<std::vector<double> > &f) const;
void setup(AtomSet& atoms);
void reset(void);
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// GlobalExtForce.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: GlobalExtForce.C,v 1.1 2010-02-20 23:13:02 fgygi Exp $
#include "GlobalExtForce.h"
#include "AtomSet.h"
#include "Atom.h"
#include "Species.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void GlobalExtForce::setup(const AtomSet& atoms)
{
// All atoms are affected, no setup needed
na_ = atoms.size();
}
////////////////////////////////////////////////////////////////////////////////
double GlobalExtForce::energy(const vector<vector<double> > &r,
vector<vector<double> > &f)
{
double sum = 0.0;
for ( int is = 0; is < f.size(); is++ )
{
const int nais = f[is].size() / 3;
for ( int ia = 0; ia < nais; ia++ )
{
f[is][3*ia+0] += force_.x / na_;
f[is][3*ia+1] += force_.y / na_;
f[is][3*ia+2] += force_.z / na_;
sum -= force_.x * r[is][3*ia+0] / na_ +
force_.y * r[is][3*ia+1] / na_ +
force_.z * r[is][3*ia+2] / na_;
}
}
return sum;
}
////////////////////////////////////////////////////////////////////////////////
ostream& GlobalExtForce::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <extforce name=\"" << name();
os << "\" type=\"" << type();
os << "\" atoms=\"_all_\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << force_.x << " "
<< setprecision(6) << force_.y << " "
<< setprecision(6) << force_.z << "\"/>";
return os;
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 The Regents of the University of California
//