Commit dafe9852 by Francois Gygi

rel1_40_0


git-svn-id: http://qboxcode.org/svn/qb/trunk@588 cba15fb0-1239-40c8-b417-11db7ca47a34
parent b7dcff3a
......@@ -3,7 +3,7 @@
// AtomsDyn.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomsDyn.h,v 1.3 2007-10-19 16:24:04 fgygi Exp $
// $Id: AtomsDyn.h,v 1.4 2008-03-26 04:57:54 fgygi Exp $
#ifndef ATOMSDYN_H
#define ATOMSDYN_H
......@@ -33,10 +33,10 @@ class AtomsDyn : public Var
}
string v = argv[1];
if ( !( v == "LOCKED" || v == "SD" || v == "SDA" || v == "MD" ) )
if ( !( v == "LOCKED" || v == "SD" || v == "SDA" || v == "CG" || v == "MD" ) )
{
if ( ui->onpe0() )
cout << " atoms_dyn must be LOCKED or SD or SDA or MD" << endl;
cout << " atoms_dyn must be LOCKED or SD or SDA or CG or MD" << endl;
return 1;
}
......
////////////////////////////////////////////////////////////////////////////////
//
// CGIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CGIonicStepper.C,v 1.1 2008-03-26 04:57:54 fgygi Exp $
#include "CGIonicStepper.h"
#include <cassert>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
bool wolfe1;
if ( first_step_ )
{
wolfe1 = true;
alpha_ = 1.0;
}
else
{
// Check if the first Wolfe condition is satisfied
// compute predicted decrease for previous step
// pred = pc_ * ( r0 - rc_ )
double pred = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
pred += pc_[is][i] * ( r0_[is][i] - rc_[is][i] );
}
}
const double sigma_wolfe = 0.1;
wolfe1 = ( e0 < ec_ - sigma_wolfe * pred );
//cout << "CGIonicStepper: pred = " << pred << endl;
//cout << "CGIonicStepper: ec: " << ec_ << endl;
//cout << "CGIonicStepper: required energy: " << ec_-sigma_wolfe * pred
// << " actual: " << e0 << endl;
//cout << "CGIonicStepper: wolfe1 = " << wolfe1 << endl;
}
if ( wolfe1 )
{
// actual decrease from last step is sufficient
// accept r0 as new current point
rc_ = r0_;
ec_ = e0;
// define new descent direction
if ( first_step_ )
{
pc_ = f0;
}
else
{
// Polak-Ribiere definition
double num = 0.0, den = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
const double fctmp = fc_[is][i];
const double f0tmp = f0[is][i];
num += f0tmp * ( f0tmp - fctmp );
den += fctmp * fctmp;
}
}
const double beta = den > 0.0 ? num/den : 0.0;
//cout << "CGIonicStepper: beta = " << beta << endl;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
pc_[is][i] = beta * pc_[is][i] + f0[is][i];
}
}
}
fc_ = f0;
alpha_ *= 1.05;
}
else
{
// backtrack
alpha_ *= 0.8;
}
//cout << "CGIonicStepper: alpha = " << alpha_ << endl;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
}
}
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
first_step_ = false;
}
////////////////////////////////////////////////////////////////////////////////
//
// CGIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CGIonicStepper.h,v 1.1 2008-03-26 04:57:54 fgygi Exp $
#ifndef CGIONICSTEPPER_H
#define CGIONICSTEPPER_H
#include "IonicStepper.h"
#include <vector>
class CGIonicStepper : public IonicStepper
{
private:
bool first_step_;
std::vector<std::vector< double> > rc_;
std::vector<std::vector< double> > pc_;
std::vector<std::vector< double> > fc_;
double ec_;
double alpha_;
public:
CGIonicStepper(Sample& s) : IonicStepper(s), first_step_(true) {}
void compute_r(double e0, const std::vector<std::vector< double> >& f0);
void compute_v(double e0, const std::vector<std::vector< double> >& f0) {}
};
#endif
......@@ -3,7 +3,7 @@
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.17 2008-03-21 00:27:11 fgygi Exp $
// $Id: MDIonicStepper.C,v 1.18 2008-03-26 04:57:54 fgygi Exp $
#include "MDIonicStepper.h"
using namespace std;
......@@ -76,7 +76,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
// if th_time_ is zero or less than dt, collision probability is one
const double collision_probability = th_time_ == 0 ? 1.0 :
const double collision_probability = th_time_ == 0 ? 1.0 :
min(fabs(dt_) / th_time_, 1.0);
if ( s_.ctxt_.onpe0() )
......@@ -128,7 +128,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
{
if ( nat > 1 )
{
collision_probability = min(1.0,fabs(dt_)/(0.5*(nat-1)*th_time_));
collision_probability = min(1.0,fabs(dt_)/(0.5*(nat-1)*th_time_));
}
}
......@@ -175,11 +175,11 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
D3vector dv12 = mu * ( lambda - v12*e12 ) * e12;
D3vector v1p = v1 + (1.0/m1) * dv12;
D3vector v2p = v2 - (1.0/m2) * dv12;
v0_[is1][3*ia1+0] = v1p.x;
v0_[is1][3*ia1+1] = v1p.y;
v0_[is1][3*ia1+2] = v1p.z;
v0_[is2][3*ia2+0] = v2p.x;
v0_[is2][3*ia2+1] = v2p.y;
v0_[is2][3*ia2+2] = v2p.z;
......
#-------------------------------------------------------------------------------
# $Id: Makefile,v 1.51 2008-03-21 00:33:09 fgygi Exp $
# $Id: Makefile,v 1.52 2008-03-26 04:57:54 fgygi Exp $
#------------------------------------------------------------------------------
#
include $(TARGET).mk
......@@ -29,7 +29,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
PSDWavefunctionStepper.o PSDAWavefunctionStepper.o \
SDCellStepper.o ConfinementPotential.o Preconditioner.o \
release.o qbox_xmlns.o isodate.o \
AndersonMixer.o SDAIonicStepper.o \
AndersonMixer.o SDAIonicStepper.o CGIonicStepper.o \
ConstraintSet.o Constraint.o DistanceConstraint.o \
AngleConstraint.o TorsionConstraint.o jacobi.o \
SampleWriter.o ComputeMLWFCmd.o BasisMapping.o MLWFTransform.o \
......@@ -226,8 +226,8 @@ BOSampleStepper.o: SlaterDet.h Basis.h Matrix.h WavefunctionStepper.h
BOSampleStepper.o: SDWavefunctionStepper.h PSDWavefunctionStepper.h
BOSampleStepper.o: PSDAWavefunctionStepper.h SDIonicStepper.h IonicStepper.h
BOSampleStepper.o: Species.h SDAIonicStepper.h AndersonMixer.h
BOSampleStepper.o: MDIonicStepper.h SDCellStepper.h CellStepper.h
BOSampleStepper.o: Preconditioner.h
BOSampleStepper.o: CGIonicStepper.h MDIonicStepper.h SDCellStepper.h
BOSampleStepper.o: CellStepper.h Preconditioner.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 Wavefunction.h
......@@ -243,6 +243,12 @@ CellMass.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
CellMass.o: ConstraintSet.h Wavefunction.h Control.h
CellStepper.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
CellStepper.o: ConstraintSet.h Wavefunction.h Control.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: Wavefunction.h Control.h Species.h
CGIonicStepper.o: IonicStepper.h Sample.h AtomSet.h Context.h Atom.h
CGIonicStepper.o: D3vector.h UnitCell.h ConstraintSet.h Wavefunction.h
CGIonicStepper.o: Control.h Species.h
ChargeDensity.o: ChargeDensity.h Timer.h Context.h Basis.h D3vector.h
ChargeDensity.o: UnitCell.h Wavefunction.h FourierTransform.h SlaterDet.h
ChargeDensity.o: Matrix.h
......
......@@ -3,7 +3,7 @@
// SDAIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.C,v 1.5 2007-10-19 16:24:04 fgygi Exp $
// $Id: SDAIonicStepper.C,v 1.6 2008-03-26 04:57:54 fgygi Exp $
#include "SDAIonicStepper.h"
#include <iostream>
......@@ -36,10 +36,19 @@ void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
<< "</sda_residual>" << endl;
mixer_.update(&f_[0],&theta_,&fbar_[0]);
if ( e0 > em_ )
{
theta_ = 0.0;
mixer_.restart();
if ( s_.ctxt_.onpe0() )
cout << " SDAIonicStepper: sda_theta = " << theta_
<< " reset to 0" << endl;
}
if ( s_.ctxt_.onpe0() )
cout << " <sda_theta> " << theta_
<< "</sda_theta>" << endl;
k = 0;
for ( int is = 0; is < r0_.size(); is++ )
{
......@@ -56,4 +65,5 @@ void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
em_ = e0;
}
......@@ -3,7 +3,7 @@
// SDAIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.h,v 1.8 2007-10-19 16:24:04 fgygi Exp $
// $Id: SDAIonicStepper.h,v 1.9 2008-03-26 04:57:54 fgygi Exp $
#ifndef SDAIONICSTEPPER_H
#define SDAIONICSTEPPER_H
......@@ -20,6 +20,7 @@ class SDAIonicStepper : public IonicStepper
double theta_;
bool first_step_;
AndersonMixer mixer_;
double em_;
public:
......@@ -29,8 +30,9 @@ class SDAIonicStepper : public IonicStepper
f_.resize(3*atoms_.size());
fbar_.resize(3*atoms_.size());
rm_ = r0_;
mixer_.set_theta_max(1.5);
mixer_.set_theta_max(1.0);
mixer_.set_theta_nc(1.0);
em_ = 1.0e50;
}
void compute_r(double e0, const std::vector<std::vector< double> >& f0);
......
......@@ -3,7 +3,7 @@
// SDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.C,v 1.4 2007-03-17 01:14:00 fgygi Exp $
// $Id: SDIonicStepper.C,v 1.5 2008-03-26 04:57:54 fgygi Exp $
#include "SDIonicStepper.h"
using namespace std;
......@@ -11,17 +11,61 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
// Steepest descent step
bool wolfe1;
if ( first_step_ )
{
wolfe1 = true;
alpha_ = 1.0;
}
else
{
// Check if the first Wolfe condition is satisfied
// compute predicted decrease for previous step
// pred = pc_ * ( r0 - rc_ )
double pred = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
pred += pc_[is][i] * ( r0_[is][i] - rc_[is][i]);
}
}
const double sigma_wolfe = 0.1;
wolfe1 = ( e0 < ec_ - sigma_wolfe * pred );
//cout << "SDIonicStepper: pred = " << pred << endl;
//cout << "SDIonicStepper: ec: " << ec_ << endl;
//cout << "SDIonicStepper: required energy: " << ec_-sigma_wolfe * pred
// << " actual: " << e0 << endl;
//cout << "SDIonicStepper: wolfe1 = " << wolfe1 << endl;
}
if ( wolfe1 )
{
// actual decrease from last step is sufficient
// accept r0 as new current point
rc_ = r0_;
ec_ = e0;
// define new descent direction
pc_ = f0;
alpha_ *= 1.05;
}
else
{
// backtrack
alpha_ *= 0.8;
}
//cout << "SDIonicStepper: alpha = " << alpha_ << endl;
for ( int is = 0; is < r0_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = r0_[is][i] + dt2bym * f0[is][i];
rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
}
}
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
first_step_ = false;
}
......@@ -3,7 +3,7 @@
// SDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.h,v 1.6 2007-10-19 16:24:04 fgygi Exp $
// $Id: SDIonicStepper.h,v 1.7 2008-03-26 04:57:54 fgygi Exp $
#ifndef SDIONICSTEPPER_H
#define SDIONICSTEPPER_H
......@@ -15,9 +15,15 @@ class SDIonicStepper : public IonicStepper
{
private:
bool first_step_;
std::vector<std::vector< double> > rc_;
std::vector<std::vector< double> > pc_;
double ec_;
double alpha_;
public:
SDIonicStepper(Sample& s) : IonicStepper(s) {}
SDIonicStepper(Sample& s) : IonicStepper(s), first_step_(true) {}
void compute_r(double e0, const std::vector<std::vector< double> >& f0);
void compute_v(double e0, const std::vector<std::vector< double> >& f0) {}
......
--------------------------------------------------------------------------------
rel1_40_0
AtomsDyn.h,CGIonicStepper.[Ch]: added CG ionic stepper.
BOSampleStepper.C: enabled niter=0 in run command.
MDIonicStepper.C: removed trailing blanks.
SDAIonicStepper: reduced theta_max, test for energy increase.
SDIonicStepper: implemented line search, removed dt, mass dependency.
--------------------------------------------------------------------------------
rel1_39_0
SampleReader.C, SpeciesReader.C: set validation scheme to Val_Auto to allow
successful parsing in the absence of a local copy of sample.xsd or species.xsd.
......
......@@ -3,10 +3,10 @@
// release.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: release.C,v 1.47 2008-03-21 00:33:09 fgygi Exp $
// $Id: release.C,v 1.48 2008-03-26 04:57:54 fgygi Exp $
#include "release.h"
std::string release(void)
{
return std::string("1.39.0");
return std::string("1.40.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