Commit 4e7b104e by Francois Gygi

added randomize_v command

git-svn-id: http://qboxcode.org/svn/qb/trunk@972 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 43035b79
......@@ -20,6 +20,7 @@
#include "AtomSet.h"
#include "Species.h"
#include "NameOf.h"
#include "sampling.h"
#include <iostream>
#include <algorithm>
#include <string>
......@@ -409,13 +410,39 @@ void AtomSet::rescale_velocities(double fac)
get_velocities(v);
for ( int is = 0; is < v.size(); is++ )
{
for ( int ia = 0; ia < v[is].size(); ia++ )
v[is][ia] *= fac;
for ( int i = 0; i < v[is].size(); i++ )
v[is][i] *= fac;
}
set_velocities(v);
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::randomize_velocities(double temp)
{
// initialize velocities with random numbers from a Maxwell-Boltzmann
// distribution at temperature temp
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
vector<vector<double> > v;
get_velocities(v);
for ( int is = 0; is < v.size(); is++ )
{
const double m = species_list[is]->mass() * 1822.89;
const double width = sqrt( boltz * temp / m );
for ( int ia = 0; ia < atom_list[is].size(); ia++ )
{
// draw pairs of unit variance gaussian random variables
double xi0, xi1, xi2, xi3; // xi3 not used
normal_dev(&xi0,&xi1);
normal_dev(&xi2,&xi3);
v[is][3*ia+0] = width * xi0;
v[is][3*ia+1] = width * xi1;
v[is][3*ia+2] = width * xi2;
}
}
set_velocities(v);
sync();
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::vcm(void) const
{
D3vector mvsum;
......
......@@ -83,6 +83,7 @@ class AtomSet
void sync(void);
void reset_velocities(void);
void rescale_velocities(double fac);
void randomize_velocities(double temp);
D3vector vcm(void) const;
D3vector dipole(void) const;
void reset_vcm(void);
......
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// RandomizeVCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef RANDOMIZEVCMD_H
#define RANDOMIZEVCMD_H
#include <iostream>
#include "UserInterface.h"
#include "Sample.h"
#include <cstdlib>
class RandomizeVCmd : public Cmd
{
public:
Sample *s;
RandomizeVCmd(Sample *sample) : s(sample) {};
char *name(void) const { return "randomize_v"; }
char *help_msg(void) const
{
return
"\n randomize_v\n\n"
" syntax: randomize_v temp\n\n"
" The randomize_v command initializes atomic velocities\n"
" with random numbers drawn from a Maxwell-Boltzmann distribution\n"
" at temperature temp\n\n";
}
int action(int argc, char **argv)
{
if ( argc != 2 )
{
if ( ui->onpe0() )
{
cout << " use: randomize_v temp" << endl;
}
return 1;
}
const double temp = atof(argv[1]);
if ( temp <= 0.0 )
{
if ( ui->onpe0() )
cout << " randomize_v: temperature must be positive" << endl;
return 1;
}
s->atoms.randomize_velocities(temp);
return 0;
}
};
#endif
......@@ -60,6 +60,7 @@ using namespace std;
#include "PlotCmd.h"
#include "PrintCmd.h"
#include "QuitCmd.h"
#include "RandomizeVCmd.h"
#include "RandomizeWfCmd.h"
#include "ResetVcmCmd.h"
#include "RescaleVCmd.h"
......@@ -253,6 +254,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new PlotCmd(s));
ui.addCmd(new PrintCmd(s));
ui.addCmd(new QuitCmd(s));
ui.addCmd(new RandomizeVCmd(s));
ui.addCmd(new RandomizeWfCmd(s));
ui.addCmd(new ResetVcmCmd(s));
ui.addCmd(new RescaleVCmd(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