Commit 9b68a022 by Francois Gygi

new IonicStepper logic.


git-svn-id: http://qboxcode.org/svn/qb/trunk@413 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 22f55599
......@@ -3,7 +3,7 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.9 2005-04-29 18:14:35 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.10 2005-06-27 22:30:28 fgygi Exp $
#include "CPSampleStepper.h"
#include "SlaterDet.h"
......@@ -11,6 +11,7 @@
#include "MDIonicStepper.h"
#include "SDCellStepper.h"
#include "Basis.h"
#include "Species.h"
#include <iostream>
#include <iomanip>
......@@ -52,6 +53,7 @@ CPSampleStepper::~CPSampleStepper(void)
////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(int niter)
{
const bool onpe0 = s_.ctxt_.onpe0();
// CP dynamics is allowed only for all doubly occupied states
// check if states are all doubly occupied
const bool wf_double_occ = (s_.wf.nel() == 2 * s_.wf.nst());
......@@ -91,6 +93,9 @@ void CPSampleStepper::step(int niter)
s_.wfv->clear();
}
if ( mdionic_stepper )
mdionic_stepper->setup_constraints();
Timer tm_iter;
tmap["charge"].start();
......@@ -109,12 +114,15 @@ void CPSampleStepper::step(int niter)
tm_iter.start();
if ( s_.ctxt_.mype() == 0 )
cout << " <iteration count=\"" << iter+1 << "\">\n";
if ( mdionic_stepper )
atoms.sync();
mdwf_stepper->update(dwf);
ekin_e = mdwf_stepper->ekin();
if ( s_.ctxt_.onpe0() )
if ( onpe0 )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
......@@ -148,13 +156,12 @@ void CPSampleStepper::step(int niter)
{
if ( iter > 0 )
{
mdionic_stepper->compute_v0(fion);
mdionic_stepper->update_v();
mdionic_stepper->compute_v(energy,fion);
}
mdionic_stepper->compute_rp(fion);
mdionic_stepper->compute_r(energy,fion);
ekin_ion = mdionic_stepper->ekin();
if ( s_.ctxt_.onpe0() )
if ( onpe0 )
{
for ( int is = 0; is < atoms.atom_list.size(); is++ )
{
......@@ -183,9 +190,23 @@ void CPSampleStepper::step(int niter)
}
}
}
#if 0
if ( s_.constraints.size() > 0 )
{
const double projected_force =
s_.constraints.projection(mdionic_stepper->r0(), fion);
if ( onpe0 )
{
cout << " <constraint_projected_force> "
<< projected_force
<< " </constraint_projected_force>"
<< endl;
}
}
#endif
}
if ( s_.ctxt_.onpe0() )
if ( onpe0 )
{
cout << " <ekin_e> " << ekin_e << " </ekin_e>\n";
......@@ -201,7 +222,7 @@ void CPSampleStepper::step(int niter)
if ( compute_stress )
{
if ( s_.ctxt_.onpe0() )
if ( onpe0 )
{
cout << "<unit_cell>" << endl;
cout << s_.wf.cell();
......@@ -224,7 +245,6 @@ void CPSampleStepper::step(int niter)
if ( compute_forces )
{
mdionic_stepper->update_r();
ef_.atoms_moved();
}
......@@ -252,6 +272,8 @@ void CPSampleStepper::step(int niter)
<< " : " << setprecision(3) << setw(9) << tmin
<< " " << setprecision(3) << setw(9) << tmax << " -->" << endl;
}
if ( compute_forces )
s_.constraints.update_constraints(dt);
} // iter
// dwf and fion now contain the forces on wavefunctions and ions at the
......@@ -262,8 +284,7 @@ void CPSampleStepper::step(int niter)
if ( compute_forces )
{
// Note: next line function call updates velocities in the AtomSet
mdionic_stepper->compute_v0(fion);
mdionic_stepper->update_v();
mdionic_stepper->compute_v(energy,fion);
}
if ( cell_stepper != 0 ) delete cell_stepper;
......
......@@ -3,35 +3,53 @@
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.8 2005-05-31 18:16:02 fgygi Exp $
// $Id: MDIonicStepper.C,v 1.9 2005-06-27 22:28:33 fgygi Exp $
#include "MDIonicStepper.h"
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_rp(const vector<vector< double> >& f0)
void MDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
// f0 contains forces at r0
// Compute new positions rp using the Verlet algorithm
// compute velocity at t+1/2 vhalf_
// Compute new positions rp using the velocity Verlet algorithm
// enforce constraints for rp
// update rm <- r0, r0 <- rp, and update atomset
ekin_ = 0.0;
atoms_.get_velocities(v0_);
for ( int is = 0; is < v0_.size(); is++ )
// compute rp
for ( int is = 0; is < r0_.size(); is++ )
{
const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
if ( dt_ != 0.0 )
const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < r0_[is].size(); i++ )
{
for ( int i = 0; i < v0_[is].size(); i++ )
{
const double v = v0_[is][i] + dtby2m * f0[is][i];
ekin_ += 0.5 * pmass_[is] * v0_[is][i] * v0_[is][i];
vhalf_[is][i] = v;
}
rp_[is][i] = r0_[is][i] + v0_[is][i] * dt_ + dt2by2m * f0[is][i];
}
}
// ekin_ is the kinetic energy computed from v0
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
{
// compute velocities v0_ using r0_, rm_ and f0(r0)
// enforce constraints for vp
// adjust velocities with the thermostat
assert(dt_ > 0.0);
for ( int is = 0; is < v0_.size(); is++ )
{
assert(pmass_[is] > 0.0);
const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < v0_[is].size(); i++ )
{
const double vhalf = ( r0_[is][i] - rm_[is][i] ) / dt_;
v0_[is][i] = vhalf + dtby2m * f0[is][i];
}
}
compute_ekin();
if ( thermostat_ )
{
eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_;
......@@ -47,53 +65,25 @@ void MDIonicStepper::compute_rp(const vector<vector< double> >& f0)
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
vhalf_[is][i] *= fac;
v0_[is][i] *= fac;
}
}
ekin_ *= fac * fac;
}
// compute rp
atoms_.get_positions(r0_);
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = r0_[is][i] + vhalf_[is][i] * dt_;
}
}
constraints_.enforce_v(r0_,v0_);
atoms_.set_velocities(v0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::compute_v0(const vector<vector< double> >& f0)
void MDIonicStepper::compute_ekin(void)
{
// compute velocities v0_ using vhalf_ and force f0(r0)
// Note: vhalf contains the velocity at t-1/2 since compute_v0 is called
// after a call to update_r()
ekin_ = 0.0;
for ( int is = 0; is < vhalf_.size(); is++ )
for ( int is = 0; is < v0_.size(); is++ )
{
assert(pmass_[is] > 0.0);
const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
for ( int i = 0; i < vhalf_[is].size(); i++ )
for ( int i = 0; i < v0_[is].size(); i++ )
{
const double v = vhalf_[is][i] + dtby2m * f0[is][i];
const double v = v0_[is][i];
ekin_ += 0.5 * pmass_[is] * v * v;
v0_[is][i] = v;
}
}
// ekin_ contains the kinetic energy computed from v0
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::update_r(void)
{
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void MDIonicStepper::update_v(void)
{
atoms_.set_velocities(v0_);
}
......@@ -3,30 +3,35 @@
// MDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.h,v 1.5 2004-12-17 23:39:06 fgygi Exp $
// $Id: MDIonicStepper.h,v 1.6 2005-06-27 22:25:34 fgygi Exp $
//
// IonicStepper is used in the following way
//
// input: r0,v0
//
// compute forces f0(r0)
// for ( k=0; k<n; k++ )
// {
// compute forces f0(r0)
// if ( k > 0 ) stepper->compute_v0(f0);
// optional: modify velocities
// stepper->update_v();
// // new r0,v0,f0 known
// stepper->compute_rp(f0); // using r0, v0 and f0
// optional: restore constraints using rp and r0
// stepper->update_r(); // r0 <- rp
// // Note: r0 and v0 are not consistent at this point
// }
// compute forces f0(r0)
// stepper->compute_v0(f0);
// stepper->update_v();
// // r0,v0,f0 known
// stepper->compute_r(f0)
// {
// computes rp using r0, v0 and f0
// restores constraints on rp using rp and r0
// updates rm<-r0, r0<-rp and update atomset positions
// }
//
// // r0,v0,f0 consistent at this point
// compute f0(r0)
//
// stepper->compute_v(f0)
// {
// computes v0 using r0,rm,f0
// restores constraints on v0 using r0, v0
// modifies velocities using the thermostat (rescaling)
// updates atomset velocities
// }
// }
// r0,v0,f0 consistent at this point
//
#ifndef MDIONICSTEPPER_H
......@@ -44,7 +49,7 @@ class MDIonicStepper : public IonicStepper
double ekin_;
double eta_;
bool thermostat_;
vector<vector< double> > vhalf_; // vhalf_[nsp_][3*na_]: v(t+dt/2)
void compute_ekin(void);
public:
......@@ -56,25 +61,19 @@ class MDIonicStepper : public IonicStepper
th_width_ = s.ctrl.th_width;
eta_ = 0.0;
ekin_ = 0.0;
vhalf_.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
{
const int nais = atoms_.na(is);
vhalf_[is].resize(3*nais);
}
atoms_.get_positions(r0_);
atoms_.get_velocities(v0_);
compute_ekin();
}
void compute_rp(const vector<vector< double> >& f0);
void compute_v0(const vector<vector< double> >& f0);
void update_r(void);
void update_v(void);
void compute_r(double e0, const vector<vector< double> >& f0);
void compute_v(double e0, const vector<vector< double> >& f0);
double eta(void) const { return eta_; }
double ekin(void) const { return ekin_; }
double temp(void) const
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
if ( ndofs_ > 0.0 )
if ( ndofs_ > 0 )
return 2.0 * ( ekin_ / boltz ) / ndofs_;
else
return 0.0;
......
......@@ -3,14 +3,14 @@
// SDAIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.C,v 1.3 2004-12-18 23:21:42 fgygi Exp $
// $Id: SDAIonicStepper.C,v 1.4 2005-06-27 22:21:32 fgygi Exp $
#include "SDAIonicStepper.h"
#include <iostream>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void SDAIonicStepper::compute_rp(const vector<vector< double> >& f0)
void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
// Steepest descent step
atoms_.get_positions(r0_);
......@@ -32,12 +32,12 @@ void SDAIonicStepper::compute_rp(const vector<vector< double> >& f0)
}
}
if ( s_.ctxt_.onpe0() )
cout << "<sda_residual> " << sum
cout << " <sda_residual> " << sum
<< "</sda_residual>" << endl;
mixer_.update(&f_[0],&theta_,&fbar_[0]);
if ( s_.ctxt_.onpe0() )
cout << "<sda_theta> " << theta_
cout << " <sda_theta> " << theta_
<< "</sda_theta>" << endl;
k = 0;
......@@ -52,18 +52,8 @@ void SDAIonicStepper::compute_rp(const vector<vector< double> >& f0)
}
}
first_step_ = false;
}
////////////////////////////////////////////////////////////////////////////////
void SDAIonicStepper::update_r(void)
{
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void SDAIonicStepper::update_v(void)
{
atoms_.reset_velocities(); // set velocities to zero
}
......@@ -3,7 +3,7 @@
// SDAIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.h,v 1.3 2004-12-18 23:21:18 fgygi Exp $
// $Id: SDAIonicStepper.h,v 1.4 2005-06-27 22:21:32 fgygi Exp $
#ifndef SDAIONICSTEPPER_H
#define SDAIONICSTEPPER_H
......@@ -15,7 +15,6 @@ class SDAIonicStepper : public IonicStepper
{
private:
vector<vector<double> > rm_;
vector<double> f_;
vector<double> fbar_;
double theta_;
......@@ -34,10 +33,8 @@ class SDAIonicStepper : public IonicStepper
//mixer_.set_theta_nc(-0.5); // default: 0.0
}
void compute_rp(const vector<vector< double> >& f0);
void compute_v0(const vector<vector< double> >& f0) {}
void update_r(void);
void update_v(void);
void compute_r(double e0, const vector<vector< double> >& f0);
void compute_v(double e0, const vector<vector< double> >& f0) {}
void reset(void) { first_step_ = true; }
};
......
......@@ -3,15 +3,14 @@
// SDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.C,v 1.2 2004-03-11 21:52:32 fgygi Exp $
// $Id: SDIonicStepper.C,v 1.3 2005-06-27 22:20:41 fgygi Exp $
#include "SDIonicStepper.h"
////////////////////////////////////////////////////////////////////////////////
void SDIonicStepper::compute_rp(const vector<vector< double> >& f0)
void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
// Steepest descent step
atoms_.get_positions(r0_);
for ( int is = 0; is < r0_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
......@@ -20,17 +19,8 @@ void SDIonicStepper::compute_rp(const vector<vector< double> >& f0)
rp_[is][i] = r0_[is][i] + dt2bym * f0[is][i];
}
}
}
////////////////////////////////////////////////////////////////////////////////
void SDIonicStepper::update_r(void)
{
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
}
////////////////////////////////////////////////////////////////////////////////
void SDIonicStepper::update_v(void)
{
atoms_.reset_velocities(); // set velocities to zero
}
......@@ -3,7 +3,7 @@
// SDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.h,v 1.3 2004-12-17 23:37:25 fgygi Exp $
// $Id: SDIonicStepper.h,v 1.4 2005-06-27 22:20:10 fgygi Exp $
#ifndef SDIONICSTEPPER_H
#define SDIONICSTEPPER_H
......@@ -18,10 +18,8 @@ class SDIonicStepper : public IonicStepper
SDIonicStepper(Sample& s) : IonicStepper(s) {}
void compute_rp(const vector<vector< double> >& f0);
void compute_v0(const vector<vector< double> >& f0) {}
void update_r(void);
void update_v(void);
void compute_r(double e0, const vector<vector< double> >& f0);
void compute_v(double e0, const vector<vector< double> >& f0) {}
};
#endif
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