Commit 2a0615be by Francois Gygi

Added CG cell stepper to optimize atomic positions and cell simultaneously.


git-svn-id: http://qboxcode.org/svn/qb/trunk@918 cba15fb0-1239-40c8-b417-11db7ca47a34
parent fc311633
......@@ -15,7 +15,6 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.58 2010-05-10 20:51:55 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -32,6 +31,7 @@
#include "MDIonicStepper.h"
#include "BMDIonicStepper.h"
#include "SDCellStepper.h"
#include "CGCellStepper.h"
#include "Preconditioner.h"
#include "AndersonMixer.h"
#include "MLWFTransform.h"
......@@ -253,6 +253,8 @@ void BOSampleStepper::step(int niter)
CellStepper* cell_stepper = 0;
if ( cell_dyn == "SD" )
cell_stepper = new SDCellStepper(s_);
else if ( cell_dyn == "CG" )
cell_stepper = new CGCellStepper(s_);
// Allocate wavefunction velocity if not available
if ( atoms_move && extrapolate_wf )
......@@ -473,7 +475,7 @@ void BOSampleStepper::step(int niter)
if ( cell_moves )
{
cell_stepper->compute_new_cell(sigma);
cell_stepper->compute_new_cell(energy,sigma,fion);
// Update cell
cell_stepper->update_cell();
......@@ -791,7 +793,7 @@ void BOSampleStepper::step(int niter)
cd_.rhog[0] = rhog_in;
cd_.update_rhor();
}
} // if nite > 1
} // if nite_ > 1
ef_.update_vhxc();
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2011 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CGCellStepper.C
//
////////////////////////////////////////////////////////////////////////////////
#include "CGCellStepper.h"
#include <valarray>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
CGCellStepper::CGCellStepper(Sample& s) : CellStepper(s), first_step_(true),
sigma1_(0.1), sigma2_(0.5)
{
nat_ = atoms_.size();
xc_.resize(3*nat_+9);
pc_.resize(3*nat_+9);
fc_.resize(3*nat_+9);
rc_.resize(atoms_.nsp());
rp_.resize(atoms_.nsp());
for ( int is = 0; is < rc_.size(); is++ )
{
rc_[is].resize(3*atoms_.na(is));
rp_[is].resize(3*atoms_.na(is));
}
linmin_.set_sigma1(sigma1_);
}
////////////////////////////////////////////////////////////////////////////////
void CGCellStepper::compute_new_cell(double e0, const valarray<double>& sigma,
const std::vector<std::vector< double> >& fion)
{
// compute new cell and ionic positions using the stress tensor sigma
// and the forces on ions fion
const UnitCell& cell = s_.wf.cell();
// total number of dofs: 3* natoms + cell parameters
valarray<double> f0(3*nat_+9), x0(3*nat_+9), xp(3*nat_+9);
// copy current positions into x0
vector<vector<double> > r0;
atoms_.get_positions(r0);
double tmp3[3];
// convert position from r0 to tau (relative) coordinates: tau = A^-1 R
int i = 0;
for ( int is = 0; is < r0.size(); is++ )
{
for ( int ia = 0; ia < r0[is].size()/3; ia++ )
{
cell.vecmult3x3(cell.amat_inv(),&r0[is][3*ia],tmp3);
x0[i++]=tmp3[0];
x0[i++]=tmp3[1];
x0[i++]=tmp3[2];
}
}
// convert forces on positions to forces on tau coordinates, store in f0
// f = A * fion
for ( int is = 0; is < fion.size(); is++ )
{
for ( int ia = 0; ia < fion[is].size()/3; ia++ )
{
cell.vecmult3x3(cell.amat(),&fion[is][3*ia],tmp3);
f0[i++]=tmp3[0];
f0[i++]=tmp3[1];
f0[i++]=tmp3[2];
}
}
// convert descent direction dcell to cell space from the stress tensor
// dcell = - omega * sigma * A
valarray<double> dcell(9); // descent direction in cell space
assert(sigma.size()==6);
// next line: local copy of sigma to circumvent compiler error
valarray<double> sigma_loc(sigma);
cell.smatmult3x3(&sigma_loc[0],cell.amat(),&dcell[0]);
// dcell *= -cell.volume();
for ( int i = 0; i < dcell.size(); i++ )
f0[3*nat_+i] = dcell[i];
// cout << "dcell = " << endl;
// cout << dcell[0] << " " << dcell[1] << " " << dcell[2] << endl;
// cout << dcell[3] << " " << dcell[4] << " " << dcell[5] << endl;
// cout << dcell[6] << " " << dcell[7] << " " << dcell[8] << endl;
// the vector f0 now contains the derivatives of the energy in tau+cell space
double fp0; // derivative f'(x) at x=xp
bool wolfe1, wolfe2; // First and second Wolfe conditions in line search
// CG algorithm
if ( !first_step_ )
{
wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
// compute fp0: projection of forces on direction pc_
fp0 = 0.0;
for ( int i = 0; i < f0.size(); i++ )
fp0 -= f0[i] * pc_[i];
wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
if ( s_.ctxt_.onpe0() )
{
cout << " CGCellStepper: fpc = " << fpc_ << endl;
cout << " CGCellStepper: fp0 = " << fp0 << endl;
cout << " CGCellStepper: ec = " << ec_ << " e0 = " << e0 << endl;
cout << " CGCellStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
<< ec_ + fpc_ * sigma1_ * alpha_ << endl;
cout << " CGCellStepper: wolfe1/wolfe2 = "
<< wolfe1 << "/" << wolfe2 << endl;
}
}
if ( first_step_ || (wolfe1 && wolfe2) )
{
// set new descent direction
// pc = f0
// define new descent direction
if ( first_step_ )
{
pc_ = f0;
}
else
{
// Polak-Ribiere definition
double num = 0.0, den = 0.0;
for ( int i = 0; i < f0.size(); i++ )
{
const double fctmp = fc_[i];
const double f0tmp = f0[i];
num += f0tmp * ( f0tmp - fctmp );
den += fctmp * fctmp;
}
double beta = den > 0.0 ? num/den : 0.0;
beta = max(beta,0.0);
if ( s_.ctxt_.onpe0() )
cout << " CGCellStepper: beta = " << beta << endl;
for ( int i = 0; i < f0.size(); i++ )
pc_[i] = beta * pc_[i] + f0[i];
}
fc_ = f0;
// fpc = d_e / d_alpha in direction pc
fpc_ = 0.0;
for ( int i = 0; i < f0.size(); i++ )
fpc_ -= fc_[i] * pc_[i];
ec_ = e0;
xc_ = x0;
fp0 = fpc_;
// reset line minimizer
linmin_.reset();
}
alpha_ = linmin_.newalpha(alpha_,e0,fp0);
if ( s_.ctxt_.onpe0() )
cout << " CGCellStepper: alpha = " << alpha_ << endl;
// update
xp = xc_ + alpha_ * pc_;
// compute atomic positions rc_[is][3*ia+j] from xc_
// compute new atomic positions rp_[is][3*ia+j] from xp
for ( int is = 0, i = 0; is < fion.size(); is++ )
{
for ( int ia = 0; ia < fion[is].size()/3; ia++ )
{
tmp3[0] = xc_[i+0];
tmp3[1] = xc_[i+1];
tmp3[2] = xc_[i+2];
cell.vecmult3x3(cell.amat(),tmp3,&rc_[is][3*ia]);
tmp3[0] = xp[i+0];
tmp3[1] = xp[i+1];
tmp3[2] = xp[i+2];
cell.vecmult3x3(cell.amat(),tmp3,&rp_[is][3*ia]);
i+=3;
}
}
// enforce constraints on atomic positions
s_.constraints.enforce_r(rc_,rp_);
// save the change of cell parameters into dcell
for ( int i = 0; i < 9; i++ )
dcell[i] = xp[3*nat_+i]-xc_[3*nat_+i];
// cellp = cell + dcell
D3vector a0p = cell.a(0) + D3vector(dcell[0],dcell[1],dcell[2]);
D3vector a1p = cell.a(1) + D3vector(dcell[3],dcell[4],dcell[5]);
D3vector a2p = cell.a(2) + D3vector(dcell[6],dcell[7],dcell[8]);
cellp = UnitCell(a0p,a1p,a2p);
// check for cell_lock constraints and modify cellp if needed
enforce_constraints(cell, cellp);
first_step_ = false;
}
////////////////////////////////////////////////////////////////////////////////
void CGCellStepper::update_cell(void)
{
const UnitCell& cell = s_.wf.cell();
s_.atoms.set_positions(rp_);
s_.atoms.set_cell(cellp);
// resize wavefunction and basis sets
s_.wf.resize(cellp,s_.wf.refcell(),s_.wf.ecut());
if ( s_.wfv != 0 )
s_.wfv->resize(cellp,s_.wf.refcell(),s_.wf.ecut());
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2011 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CGCellStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef CGCELLSTEPPER_H
#define CGCELLSTEPPER_H
#include "CellStepper.h"
#include "LineMinimizer.h"
class CGCellStepper : public CellStepper
{
private:
bool first_step_;
int nat_;
std::valarray<double> xc_, pc_, fc_;
std::vector<std::vector<double> > rc_, rp_;
double ec_, fpc_;
double alpha_, sigma1_, sigma2_;
LineMinimizer linmin_;
public:
CGCellStepper(Sample& s);
void compute_new_cell(double e0, const std::valarray<double>& sigma,
const std::vector<std::vector< double> >& fion);
void update_cell(void);
double ekin(void) const { return 0.0; }
};
#endif
......@@ -15,13 +15,13 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.23 2010-02-20 23:13:02 fgygi Exp $
#include "CPSampleStepper.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
#include "SDCellStepper.h"
#include "CGCellStepper.h"
#include "Basis.h"
#include "Species.h"
......@@ -114,6 +114,8 @@ void CPSampleStepper::step(int niter)
CellStepper* cell_stepper = 0;
if ( cell_dyn == "SD" )
cell_stepper = new SDCellStepper(s_);
else if ( cell_dyn == "CG" )
cell_stepper = new CGCellStepper(s_);
if ( s_.wfv == 0 )
{
......@@ -256,7 +258,7 @@ void CPSampleStepper::step(int niter)
if ( cell_dyn != "LOCKED" )
{
cell_stepper->compute_new_cell(sigma);
cell_stepper->compute_new_cell(energy,sigma,fion);
// Update cell
cell_stepper->update_cell();
......
......@@ -15,7 +15,6 @@
// CellDyn.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CellDyn.h,v 1.5 2008-09-08 15:56:18 fgygi Exp $
#ifndef CELLDYN_H
#define CELLDYN_H
......@@ -47,10 +46,10 @@ class CellDyn : public Var
}
string v = argv[1];
if ( !( v == "LOCKED" || v == "SD" ) )
if ( !( v == "LOCKED" || v == "SD" || v == "CG" ) )
{
if ( ui->onpe0() )
cout << " cell_dyn must be in [LOCKED,SD]" << endl;
cout << " cell_dyn must be LOCKED, SD or CG" << endl;
return 1;
}
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2011 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CellStepper.C
//
////////////////////////////////////////////////////////////////////////////////
#include "CellStepper.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void CellStepper::enforce_constraints(const UnitCell& cell, UnitCell& cellp)
{
// modify cellp so as to satisfy the constraints expressed in cell_lock
string cell_lock = s_.ctrl.cell_lock;
valarray<double> dcell(9); // cellp - cell
for ( int i = 0; i < 9; i++ )
dcell[i] = cellp.amat(i) - cell.amat(i);
if ( cell_lock != "OFF" )
{
// constraints on the cell derivatives
if ( cell_lock.find("A") != string::npos )
{
// vector A is locked
dcell[0] = dcell[1] = dcell[2] = 0.0;
}
if ( cell_lock.find("B") != string::npos )
{
// vector B is locked
dcell[3] = dcell[4] = dcell[5] = 0.0;
}
if ( cell_lock.find("C") != string::npos )
{
// vector C is locked
dcell[6] = dcell[7] = dcell[8] = 0.0;
}
// Check if cell shape should be preserved (if "S" is present in cell_lock)
// The only changes allowed are renormalizations of a,b,c
if ( cell_lock.find("S") != string::npos )
{
// projection of d in the direction of e
D3vector d,e;
d = D3vector(dcell[0],dcell[1],dcell[2]);
e = cell.a(0) / length(cell.a(0));
d = (d * e) * e;
dcell[0] = d.x; dcell[1] = d.y; dcell[2] = d.z;
d = D3vector(dcell[3],dcell[4],dcell[5]);
e = cell.a(1) / length(cell.a(1));
d = (d * e) * e;
dcell[3] = d.x; dcell[4] = d.y; dcell[5] = d.z;
d = D3vector(dcell[6],dcell[7],dcell[8]);
e = cell.a(2) / length(cell.a(2));
d = (d * e) * e;
dcell[6] = d.x; dcell[7] = d.y; dcell[8] = d.z;
}
if ( cell_lock == "R" )
{
// preserve aspect ratio
// dcell must be proportional to A
// All vectors are rescaled by the same constant
// rescale cell by coefficient alpha, i.e.
// dcell = alpha * A
// where alpha * A is the projection of dcell in the direction of A
// alpha = tr (A^T * dcell) / || A ||^2 (matrix scalar product)
const double *a = cell.amat();
const double num = a[0]*dcell[0] + a[1]*dcell[1] + a[2]*dcell[2] +
a[3]*dcell[3] + a[4]*dcell[4] + a[5]*dcell[5] +
a[6]*dcell[6] + a[7]*dcell[7] + a[8]*dcell[8];
const double denom = a[0]*a[0] + a[1]*a[1] + a[2]*a[2] +
a[3]*a[3] + a[4]*a[4] + a[5]*a[5] +
a[6]*a[6] + a[7]*a[7] + a[8]*a[8];
const double alpha = num / denom;
dcell = valarray<double>(a,9);
dcell *= alpha;
}
}
// dcell now represents a change compatible with constraints
// cellp = cell + dcell
D3vector a0p = cell.a(0) + D3vector(dcell[0],dcell[1],dcell[2]);
D3vector a1p = cell.a(1) + D3vector(dcell[3],dcell[4],dcell[5]);
D3vector a2p = cell.a(2) + D3vector(dcell[6],dcell[7],dcell[8]);
cellp = UnitCell(a0p,a1p,a2p);
}
......@@ -15,7 +15,6 @@
// CellStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CellStepper.h,v 1.5 2008-09-08 15:56:18 fgygi Exp $
#ifndef CELLSTEPPER_H
#define CELLSTEPPER_H
......@@ -36,7 +35,9 @@ class CellStepper
CellStepper (Sample& s) : s_(s), atoms_(s.atoms), ekin_(0.0) {}
virtual void compute_new_cell(const std::valarray<double>& sigma) = 0;
virtual void compute_new_cell(double e0,const std::valarray<double>& sigma,
const std::vector<std::vector< double> >& f0) = 0;
void enforce_constraints(const UnitCell& cell, UnitCell& cellp);
virtual void update_cell(void) = 0;
double ekin(void) const { return ekin_; }
......
......@@ -40,8 +40,8 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
MDWavefunctionStepper.o SDIonicStepper.o MDIonicStepper.o \
BMDIonicStepper.o \
PSDWavefunctionStepper.o PSDAWavefunctionStepper.o \
JDWavefunctionStepper.o \
SDCellStepper.o ConfinementPotential.o Preconditioner.o \
JDWavefunctionStepper.o CellStepper.o SDCellStepper.o \
CGCellStepper.o ConfinementPotential.o Preconditioner.o \
release.o qbox_xmlns.o isodate.o \
AndersonMixer.o SDAIonicStepper.o CGIonicStepper.o \
ConstraintSet.o Constraint.o PositionConstraint.o DistanceConstraint.o \
......@@ -211,7 +211,7 @@ SampleStepper.o EnergyFunctional.o
rm -f libqb.a
#------------------------------------------------------------------------------
ctags :
etags -o tags *.[Ch]
ctags -o tags *.[Ch]
#------------------------------------------------------------------------------
html :
enscript -Ecpp --color -Whtml --toc -pqbsrc.html *.h *.C
......@@ -269,8 +269,9 @@ BOSampleStepper.o: PSDWavefunctionStepper.h PSDAWavefunctionStepper.h
BOSampleStepper.o: JDWavefunctionStepper.h SDIonicStepper.h IonicStepper.h
BOSampleStepper.o: Species.h SDAIonicStepper.h LineMinimizer.h
BOSampleStepper.o: CGIonicStepper.h MDIonicStepper.h BMDIonicStepper.h
BOSampleStepper.o: SDCellStepper.h CellStepper.h Preconditioner.h
BOSampleStepper.o: AndersonMixer.h MLWFTransform.h BasisMapping.h
BOSampleStepper.o: SDCellStepper.h CellStepper.h CGCellStepper.h
BOSampleStepper.o: Preconditioner.h AndersonMixer.h MLWFTransform.h
BOSampleStepper.o: BasisMapping.h
BOSampleStepper.o: SampleStepper.h Timer.h EnergyFunctional.h ChargeDensity.h
BOSampleStepper.o: Context.h StructureFactor.h Sample.h AtomSet.h Atom.h
BOSampleStepper.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
......@@ -284,8 +285,17 @@ CellLock.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
CellLock.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
CellMass.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
CellMass.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
CellStepper.o: CellStepper.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
CellStepper.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h
CellStepper.o: Control.h
CellStepper.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
CellStepper.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
CGCellStepper.o: CGCellStepper.h CellStepper.h Sample.h AtomSet.h Context.h
CGCellStepper.o: Atom.h D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
CGCellStepper.o: Wavefunction.h Control.h LineMinimizer.h
CGCellStepper.o: CellStepper.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
CGCellStepper.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h
CGCellStepper.o: Control.h LineMinimizer.h
CGIonicStepper.o: CGIonicStepper.h IonicStepper.h Sample.h AtomSet.h
CGIonicStepper.o: Context.h Atom.h D3vector.h UnitCell.h ConstraintSet.h
CGIonicStepper.o: ExtForceSet.h Wavefunction.h Control.h Species.h
......@@ -327,7 +337,7 @@ CPSampleStepper.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h
CPSampleStepper.o: Control.h SlaterDet.h Basis.h Matrix.h
CPSampleStepper.o: MDWavefunctionStepper.h WavefunctionStepper.h
CPSampleStepper.o: MDIonicStepper.h IonicStepper.h Species.h SDCellStepper.h
CPSampleStepper.o: CellStepper.h
CPSampleStepper.o: CellStepper.h CGCellStepper.h LineMinimizer.h
CPSampleStepper.o: SampleStepper.h Timer.h EnergyFunctional.h ChargeDensity.h
CPSampleStepper.o: Context.h StructureFactor.h Sample.h AtomSet.h Atom.h
CPSampleStepper.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
......
......@@ -15,15 +15,16 @@
// SDCellStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDCellStepper.C,v 1.11 2008-09-26 21:05:18 fgygi Exp $
#include "SDCellStepper.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void SDCellStepper::compute_new_cell(const valarray<double>& sigma)
void SDCellStepper::compute_new_cell(double e0, const valarray<double>& sigma,
const std::vector<std::vector< double> >& f0)
{
// compute new cell by adding correction dcell_ij
// compute new cell by adding correction dcell_ij. Steepest descent algorithm.
// forces f0 are not used
valarray<double> dcell(9);
const UnitCell& cell = s_.wf.cell();
......@@ -48,80 +49,16 @@ void SDCellStepper::compute_new_cell(const valarray<double>& sigma)
// cell.smatmult3x3(&sigma[0],cell.amat(),&dcell[0]);
dcell *= -cell.volume();
string cell_lock = s_.ctrl.cell_lock;
if ( cell_lock != "OFF" )
{
// constraints on the cell derivatives
if ( cell_lock.find("A") != string::npos )
{
// vector A is locked
dcell[0] = dcell[1] = dcell[2] = 0.0;
}
if ( cell_lock.find("B") != string::npos )
{
// vector B is locked
dcell[3] = dcell[4] = dcell[5] = 0.0;
}
if ( cell_lock.find("C") != string::npos )
{
// vector C is locked
dcell[6] = dcell[7] = dcell[8] = 0.0;
}
// Check if cell shape should be preserved (if "S" is present in cell_lock)
// The only changes allowed are renormalizations of a,b,c
if ( cell_lock.find("S") != string::npos )
{
// projection of d in the direction of e
D3vector d,e;
d = D3vector(dcell[0],dcell[1],dcell[2]);
e = cell.a(0) / length(cell.a(0));
d = (d * e) * e;
dcell[0] = d.x; dcell[1] = d.y; dcell[2] = d.z;
d = D3vector(dcell[3],dcell[4],dcell[5]);
e = cell.a(1) / length(cell.a(1));
d = (d * e) * e;
dcell[3] = d.x; dcell[4] = d.y; dcell[5] = d.z;
d = D3vector(dcell[6],dcell[7],dcell[8]);
e = cell.a(2) / length(cell.a(2));
d = (d * e) * e;
dcell[6] = d.x; dcell[7] = d.y; dcell[8] = d.z;
}
if ( cell_lock == "R" )
{
// preserve aspect ratio
// dcell must be proportional to A
// All vectors are rescaled by the same constant
// rescale cell by coefficient alpha, i.e.
// dcell = alpha * A
// where alpha * A is the projection of dcell in the direction of A
// alpha = tr (A^T * dcell) / || A ||^2 (matrix scalar product)
const double *a = cell.amat();
const double num = a[0]*dcell[0] + a[1]*dcell[1] + a[2]*dcell[2] +
a[3]*dcell[3] + a[4]*dcell[4] + a[5]*dcell[5] +
a[6]*dcell[6] + a[7]*dcell[7] + a[8]*dcell[8];
const double denom = a[0]*a[0] + a[1]*a[1] + a[2]*a[2] +
a[3]*a[3] + a[4]*a[4] + a[5]*a[5] +
a[6]*a[6] + a[7]*a[7] + a[8]*a[8];
const double alpha = num / denom;
dcell = valarray<double>(a,9);
dcell *= alpha;
}
}
// new cell: cellp = cell - dcell * dt^2 / cell_mass
const double dt = s_.ctrl.dt;
const double dt2bym = dt*dt/cell_mass;
// cellp = cell - dcell * dt^2 / cell_mass
D3vector a0p = cell.a(0) - dt2bym * D3vector(dcell[0],dcell[1],dcell[2]);
D3vector a1p = cell.a(1) - dt2bym * D3vector(dcell[3],dcell[4],dcell[5]);
D3vector a2p = cell.a(2) - dt2bym * D3vector(dcell[6],dcell[7],dcell[8]);
cellp = UnitCell(a0p,a1p,a2p);
// check for cell_lock constraints and modify cellp if needed
enforce_constraints(cell, cellp);
}
////////////////////////////////////////////////////////////////////////////////
......@@ -152,8 +89,6 @@ void SDCellStepper::update_cell(void)
s_.atoms.set_cell(cellp);
// resize wavefunction and basis sets
//cout << " SDCellStepper::update_cell" << endl;
s_.wf.resize(cellp,s_.wf.refcell(),s_.wf.ecut());
if ( s_.wfv != 0 )
s_.wfv->resize(cellp,s_.wf.refcell(),s_.wf.ecut());
......
......@@ -15,7 +15,6 @@
// SDCellStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDCellStepper.h,v 1.6 2008-09-15 14:58:52 fgygi Exp $
#ifndef SDCELLSTEPPER_H
#define SDCELLSTEPPER_H
......@@ -30,7 +29,8 @@ class SDCellStepper : public CellStepper
SDCellStepper(Sample& s) : CellStepper(s) {}
void compute_new_cell(const std::valarray<double>& sigma);
void compute_new_cell(double e0, const std::valarray<double>& sigma,
const std::vector<std::vector< double> >& f0);
void update_cell(void);
double ekin(void) const { return 0.0; }
};
......
TODO:
LineMinimizer.h: avoid getting blocked in bracketing mode
indefinitely when the accuracy of forces is insufficient.
Sample.h: replace ExtForceSet and ConstraintSet members by pointers.
Check behavior of net_charge during save/load.
--------------------------------------------------------------------------------
Issues:
BOSampleStepper.C:Using run 0 modifies the wavefunctions. Need to set
wf_dyn=LOCKED to only compute energy without modif of wfs.
kpgen for non-rectangular cells yields non-symmetric kpoint sets.
--------------------------------------------------------------------------------
rel1_54_0 candidate
LineMinimizer.h: limited the number of iterations in bracketing mode to 4
UnitCell.C: accept unit cells with orientation leading to negative volume and
define the volume as fabs(a0*(a1^a2))
CellStepper.C, CGCellStepper.C: implemented CG algorithm to optimize both
cell parameters and atomic positions. Moved implementation of common CellStepper
member functions to CellStepper.C.
--------------------------------------------------------------------------------
rel1_53_0
RseedCmd.h: rseed command to initialize random number generator
BOSampleStepper.C: reintroduced wf extrapolation when nite>1 (leads to
......
......@@ -15,10 +15,9 @@
// release.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: release.C,v 1.84 2010-05-10 20:08:12 fgygi Exp $
#include "release.h"
std::string release(void)
{
return std::string("1.53.0");
return std::string("1.54.0");
}
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