Commit 7c62abc5 by Francois Gygi

Implement reset_rotation command

parent 0b1e4fad
......@@ -503,6 +503,25 @@ void AtomSet::randomize_velocities(double temp)
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::rcm(void) const
{
D3vector mrsum;
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 r = atom_list[is][ia]->position();
mrsum += mass * r;
msum += mass;
}
}
if ( msum == 0.0 ) return D3vector(0,0,0);
return mrsum / msum;
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::vcm(void) const
{
D3vector mvsum;
......@@ -542,6 +561,142 @@ void AtomSet::reset_vcm(void)
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::reset_rotation(void)
{
// check for special case of zero or one atom
if ( size() < 2 ) return;
D3vector rc = rcm();
D3vector vc = vcm();
vector<vector<double> > rt;
get_positions(rt);
vector<vector<double> > vt;
get_velocities(vt);
// compute angular momentum w.r.t. the center of mass
D3vector L;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
D3vector v(vt[is][3*ia+0],vt[is][3*ia+1],vt[is][3*ia+2]);
v -= vc;
L += mass * ( r ^ v );
}
}
// compute inertia tensor a and vector omega
// check for special case of all atoms aligned, for which the
// inertia tensor has a zero eigenvalue
// check if all atoms are aligned
// collect all positions in a single vector of D3vectors
vector<D3vector> rv(size());
int iat = 0;
for ( int is = 0; is < vt.size(); is++ )
for ( int ia = 0; ia < na(is); ia++ )
rv[iat++] = D3vector(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
// normalized direction e = (rv[1]-rv[0])
D3vector e = normalized(rv[1]-rv[0]);
bool aligned = true;
for ( int i = 2; (i < size()) && aligned; i++ )
{
D3vector u = normalized(rv[i]-rv[0]);
aligned &= length(u^e) < 1.e-6;
}
D3vector omega;
// compute the inertia tensor
// treat separately the case of all atoms aligned
if ( aligned )
{
// inertia tensor reduces to a scalar
double a = 0.0;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
a += mass * norm2(r);
}
}
omega = L / a;
}
else
{
double a00,a01,a02,a10,a11,a12,a20,a21,a22;
a00=a01=a02=a10=a11=a12=a20=a21=a22=0.0;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
a00 += mass * ( r.y * r.y + r.z * r.z );
a11 += mass * ( r.x * r.x + r.z * r.z );
a22 += mass * ( r.x * r.x + r.y * r.y );
a01 -= mass * r.x * r.y;
a02 -= mass * r.x * r.z;
a12 -= mass * r.y * r.z;
}
}
a10 = a01;
a20 = a02;
a21 = a12;
// inverse b of the inertia tensor a
// determinant
double det = a00 * ( a11 * a22 - a21 * a12 ) -
a01 * ( a10 * a22 - a20 * a12 ) +
a02 * ( a10 * a21 - a20 * a11 );
// the determinant must be positive
assert(det>1.e-8);
double b00,b01,b02,b10,b11,b12,b20,b21,b22;
b00 = ( -a12*a21 + a11*a22 ) / det;
b10 = ( a12*a20 - a10*a22 ) / det;
b20 = ( -a11*a20 + a10*a21 ) / det;
b01 = ( a02*a21 - a01*a22 ) / det;
b11 = ( -a02*a20 + a00*a22 ) / det;
b21 = ( a01*a20 - a00*a21 ) / det;
b02 = ( -a02*a11 + a01*a12 ) / det;
b12 = ( a02*a10 - a00*a12 ) / det;
b22 = ( -a01*a10 + a00*a11 ) / det;
// omega = inverse(a) * L
omega.x = b00 * L.x + b01 * L.y + b02 * L.z;
omega.y = b10 * L.x + b11 * L.y + b12 * L.z;
omega.z = b20 * L.x + b21 * L.y + b22 * L.z;
}
// correct velocities: v = v - omega ^ r
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
D3vector dv = omega ^ r;
vt[is][3*ia+0] -= dv.x;
vt[is][3*ia+1] -= dv.y;
vt[is][3*ia+2] -= dv.z;
}
}
set_velocities(vt);
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::fold_in_ws(void)
{
vector<vector<double> > p;
......
......@@ -85,10 +85,12 @@ class AtomSet
void rescale_velocities(double fac);
void randomize_velocities(double temp);
void randomize_positions(double amplitude);
D3vector rcm(void) const;
D3vector vcm(void) const;
D3vector dipole(void) const;
D3tensor quadrupole(void) const;
void reset_vcm(void);
void reset_rotation(void);
void fold_in_ws(void);
int size(void) const;
};
......
......@@ -291,7 +291,7 @@ CellStepper.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
CellStepper.o: Wavefunction.h Control.h
ChargeDensity.o: ChargeDensity.h Timer.h Context.h blacs.h Basis.h D3vector.h
ChargeDensity.o: UnitCell.h Wavefunction.h FourierTransform.h SlaterDet.h
ChargeDensity.o: Matrix.h
ChargeDensity.o: Matrix.h blas.h
ChargeDensity.o: Timer.h Context.h blacs.h
ChargeMixCoeff.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
ChargeMixCoeff.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
......@@ -378,11 +378,13 @@ ExchangeOperator.o: VectorLess.h ExchangeOperator.h Sample.h AtomSet.h
ExchangeOperator.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h D3tensor.h
ExchangeOperator.o: blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
ExchangeOperator.o: Control.h SlaterDet.h Basis.h Matrix.h Timer.h
ExchangeOperator.o: FourierTransform.h Bisection.h
ExchangeOperator.o: FourierTransform.h InteractionPotential.h Bisection.h
ExchangeOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
ExchangeOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h
ExchangeOperator.o: ExtForceSet.h Wavefunction.h Control.h SlaterDet.h
ExchangeOperator.o: Basis.h Matrix.h Timer.h FourierTransform.h
ExchangeOperator.o: InteractionPotential.h
ExponentialIntegral.o: ExponentialIntegral.h
ExtForce.o: ExtForce.h D3vector.h
ExtForce.o: D3vector.h
ExtForceCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
......@@ -406,6 +408,9 @@ FourierTransform.o: Timer.h
GlobalExtForce.o: GlobalExtForce.h ExtForce.h D3vector.h AtomSet.h Context.h
GlobalExtForce.o: blacs.h Atom.h UnitCell.h D3tensor.h blas.h Species.h
GlobalExtForce.o: ExtForce.h D3vector.h
HSEFunctional.o: HSEFunctional.h XCFunctional.h InteractionPotential.h
HSEFunctional.o: ExponentialIntegral.h
HSEFunctional.o: XCFunctional.h InteractionPotential.h
HelpCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
HelpCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
HelpCmd.o: ExtForceSet.h Wavefunction.h Control.h
......@@ -549,6 +554,9 @@ RefCell.o: Control.h
RescaleVCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
RescaleVCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
RescaleVCmd.o: ExtForceSet.h Wavefunction.h Control.h
ResetRotationCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h
ResetRotationCmd.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
ResetRotationCmd.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
ResetVcmCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
ResetVcmCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
ResetVcmCmd.o: ExtForceSet.h Wavefunction.h Control.h
......@@ -702,14 +710,16 @@ XCOperator.o: XCOperator.h Sample.h AtomSet.h Context.h blacs.h Atom.h
XCOperator.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
XCOperator.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
XCOperator.o: XCPotential.h ExchangeOperator.h SlaterDet.h Basis.h Matrix.h
XCOperator.o: FourierTransform.h
XCOperator.o: FourierTransform.h InteractionPotential.h HSEFunctional.h
XCOperator.o: XCFunctional.h
XCOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
XCOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
XCOperator.o: Wavefunction.h Control.h
XCPotential.o: XCPotential.h Control.h D3vector.h ChargeDensity.h Timer.h
XCPotential.o: Context.h blacs.h LDAFunctional.h XCFunctional.h
XCPotential.o: VWNFunctional.h PBEFunctional.h BLYPFunctional.h
XCPotential.o: B3LYPFunctional.h Basis.h UnitCell.h FourierTransform.h blas.h
XCPotential.o: HSEFunctional.h InteractionPotential.h B3LYPFunctional.h
XCPotential.o: Basis.h UnitCell.h FourierTransform.h blas.h
XCPotential.o: Control.h D3vector.h ChargeDensity.h Timer.h Context.h blacs.h
XMLGFPreprocessor.o: Timer.h Context.h blacs.h Base64Transcoder.h Matrix.h
XMLGFPreprocessor.o: XMLGFPreprocessor.h
......@@ -728,16 +738,16 @@ qb.o: Control.h Timer.h AngleCmd.h AtomCmd.h ComputeMLWFCmd.h MLWFTransform.h
qb.o: BasisMapping.h ConstraintCmd.h DistanceCmd.h ExtForceCmd.h
qb.o: FoldInWsCmd.h HelpCmd.h KpointCmd.h ListAtomsCmd.h ListSpeciesCmd.h
qb.o: LoadCmd.h MoveCmd.h PartialChargeCmd.h PlotCmd.h PrintCmd.h QuitCmd.h
qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetVcmCmd.h
qb.o: RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h SetVelocityCmd.h
qb.o: SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h BisectionCmd.h
qb.o: Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AtomsDyn.h BlHF.h
qb.o: BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h ChargeMixCoeff.h
qb.o: ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h Ecutprec.h
qb.o: Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h FermiTemp.h
qb.o: IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h Nrowmax.h Nspin.h
qb.o: RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h
qb.o: WfDiag.h WfDyn.h Xc.h
qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetRotationCmd.h
qb.o: ResetVcmCmd.h RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h
qb.o: SetVelocityCmd.h SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h
qb.o: BisectionCmd.h Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h
qb.o: AtomsDyn.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h
qb.o: ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h
qb.o: Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h
qb.o: FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h
qb.o: Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h
qb.o: ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.h
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ResetRotationCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef RESETROTATIONCMD_H
#define RESETROTATIONCMD_H
#include <string>
#include <cstdlib>
#include <iostream>
using namespace std;
#include "UserInterface.h"
#include "Sample.h"
class ResetRotationCmd : public Cmd
{
public:
Sample *s;
ResetRotationCmd(Sample *sample) : s(sample) {};
const char *name(void) const { return "reset_rotation"; }
const char *help_msg(void) const
{
return
"\n reset_rotation\n\n"
" syntax: reset_rotation\n\n"
" The reset_rotation command adjusts velocities of the atoms\n"
" to cancel the global rotation of the system.\n\n";
}
int action(int argc, char **argv)
{
s->atoms.reset_rotation();
return 0;
}
};
#endif
......@@ -63,6 +63,7 @@ using namespace std;
#include "RandomizeRCmd.h"
#include "RandomizeVCmd.h"
#include "RandomizeWfCmd.h"
#include "ResetRotationCmd.h"
#include "ResetVcmCmd.h"
#include "RescaleVCmd.h"
#include "RseedCmd.h"
......@@ -280,6 +281,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new RandomizeVCmd(s));
ui.addCmd(new RandomizeWfCmd(s));
ui.addCmd(new RescaleVCmd(s));
ui.addCmd(new ResetRotationCmd(s));
ui.addCmd(new ResetVcmCmd(s));
ui.addCmd(new RseedCmd(s));
ui.addCmd(new RunCmd(s));
......
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