Commit 40e81ab3 by Francois Gygi

Merge branch 'develop'

parents 612ee252 beb8bee0
......@@ -56,6 +56,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
ElectricEnthalpy.o PartialChargeCmd.o \
ExternalPotential.o ResponseCmd.o \
Function3d.o Function3dHandler.o \
SpectrumCmd.o \
$(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
# to include VERSION info in release string, use:
......@@ -491,17 +492,6 @@ JDWavefunctionStepper.o: D3vector.h UnitCell.h
KpointCmd.o: UserInterface.h D3vector.h Sample.h AtomSet.h Context.h blacs.h
KpointCmd.o: Atom.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
KpointCmd.o: ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: LBFGS_IonicStepper.h IonicStepper.h Sample.h AtomSet.h
LBFGS_IonicStepper.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h
LBFGS_IonicStepper.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_IonicStepper.o: Wavefunction.h Control.h Species.h
LBFGS_IonicStepper.o: IonicStepper.h Sample.h AtomSet.h Context.h blacs.h
LBFGS_IonicStepper.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
LBFGS_IonicStepper.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: Species.h
LBFGS_hessian.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
LBFGS_hessian.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_hessian.o: Wavefunction.h Control.h
LDAFunctional.o: LDAFunctional.h XCFunctional.h
LDAFunctional.o: XCFunctional.h
LineMinimizer.o: LineMinimizer.h
......@@ -738,6 +728,13 @@ SpeciesHandler.o: SpeciesHandler.h StructureHandler.h Species.h StrX.h
SpeciesHandler.o: StructureHandler.h
SpeciesReader.o: Species.h SpeciesReader.h StructuredDocumentHandler.h StrX.h
SpeciesReader.o: StructureHandler.h SpeciesHandler.h
SpectrumCmd.o: SpectrumCmd.h UserInterface.h Context.h blacs.h
SpectrumCmd.o: ChargeDensity.h Timer.h D3vector.h EnergyFunctional.h
SpectrumCmd.o: StructureFactor.h ElectricEnthalpy.h Matrix.h Wavefunction.h
SpectrumCmd.o: UnitCell.h SlaterDet.h Basis.h Sample.h AtomSet.h Atom.h
SpectrumCmd.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Control.h
SpectrumCmd.o: MLWFTransform.h BasisMapping.h
SpectrumCmd.o: UserInterface.h
StatusCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
StatusCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
StatusCmd.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
......@@ -806,16 +803,14 @@ 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 HSEFunctional.h XCFunctional.h
XCOperator.o: RSHFunctional.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: SCANFunctional.h HSEFunctional.h RSHFunctional.h
XCPotential.o: B3LYPFunctional.h BHandHLYPFunctional.h Sample.h AtomSet.h
XCPotential.o: Atom.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
XCPotential.o: ExtForceSet.h Wavefunction.h SlaterDet.h Basis.h Matrix.h
XCPotential.o: StructureFactor.h XCOperator.h NonLocalPotential.h
XCPotential.o: ConfinementPotential.h ElectricEnthalpy.h FourierTransform.h
XCPotential.o: Control.h D3vector.h ChargeDensity.h Timer.h Context.h blacs.h
XCPotential.o: XCPotential.h ChargeDensity.h Timer.h Context.h blacs.h
XCPotential.o: D3vector.h LDAFunctional.h XCFunctional.h VWNFunctional.h
XCPotential.o: PBEFunctional.h BLYPFunctional.h HSEFunctional.h
XCPotential.o: RSHFunctional.h B3LYPFunctional.h BHandHLYPFunctional.h
XCPotential.o: SCANFunctional.h Sample.h AtomSet.h Atom.h UnitCell.h
XCPotential.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
XCPotential.o: Control.h SlaterDet.h Basis.h Matrix.h FourierTransform.h
XCPotential.o: ChargeDensity.h Timer.h Context.h blacs.h D3vector.h
XMLGFPreprocessor.o: Timer.h Context.h blacs.h Base64Transcoder.h Matrix.h
XMLGFPreprocessor.o: XMLGFPreprocessor.h
XMLGFPreprocessor.o: Matrix.h Context.h blacs.h
......@@ -834,16 +829,16 @@ 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 ResetRotationCmd.h
qb.o: ResetVcmCmd.h RescaleVCmd.h ResponseCmd.h RseedCmd.h RunCmd.h SaveCmd.h
qb.o: SetCmd.h SetVelocityCmd.h SpeciesCmd.h StatusCmd.h ChargeDensity.h
qb.o: FourierTransform.h StrainCmd.h TorsionCmd.h BisectionCmd.h Bisection.h
qb.o: SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AlphaRSH.h AtomsDyn.h
qb.o: BetaRSH.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 ForceTol.h Polarization.h Emass.h
qb.o: ExtStress.h FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h
qb.o: NetCharge.h Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h StressTol.h
qb.o: Thermostat.h ThTemp.h ThTime.h ThWidth.h Vext.h ExternalPotential.h
qb.o: WfDiag.h WfDyn.h Xc.h
qb.o: SetCmd.h SetVelocityCmd.h SpeciesCmd.h SpectrumCmd.h StatusCmd.h
qb.o: ChargeDensity.h FourierTransform.h StrainCmd.h TorsionCmd.h
qb.o: BisectionCmd.h Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h
qb.o: AlphaRSH.h AtomsDyn.h BetaRSH.h BlHF.h BtHF.h Cell.h CellDyn.h
qb.o: CellLock.h CellMass.h ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h
qb.o: Debug.h Dspin.h Ecut.h Ecutprec.h Ecuts.h Efield.h ForceTol.h
qb.o: Polarization.h Emass.h ExtStress.h FermiTemp.h IterCmd.h
qb.o: IterCmdPeriod.h Dt.h MuRSH.h Nempty.h NetCharge.h Nrowmax.h Nspin.h
qb.o: RefCell.h ScfTol.h Stress.h StressTol.h Thermostat.h ThTemp.h ThTime.h
qb.o: ThWidth.h Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.h
......@@ -887,7 +882,6 @@ testWavefunction.o: SlaterDet.h Basis.h Matrix.h Timer.h
testXCFunctional.o: LDAFunctional.h XCFunctional.h PBEFunctional.h Timer.h
testXMLGFPreprocessor.o: Context.h blacs.h Matrix.h XMLGFPreprocessor.h
test_fftw.o: Timer.h readTSC.h
test_sym.o: Basis.h D3vector.h UnitCell.h
test_vext.o: Function3d.h D3vector.h
testjacobi.o: Timer.h Context.h blacs.h Matrix.h jacobi.h
testjade.o: Timer.h Context.h blacs.h Matrix.h jade.h
......
......@@ -129,7 +129,7 @@ void SlaterDet::resize(const UnitCell& cell, const UnitCell& refcell,
// perform normal resize operations, possibly resetting contents of c_
basis_->resize(cell,refcell,ecut);
occ_.resize(nst);
eig_.resize(nst);
eig_.resize(nst,0.0);
const int mb = basis_->maxlocalsize();
const int m = ctxt_.nprow() * mb;
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2020 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// SpectrumCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
#include "SpectrumCmd.h"
#include <iostream>
#include <fstream>
#include <cassert>
#include "Context.h"
#include "ChargeDensity.h"
#include "EnergyFunctional.h"
#include "MLWFTransform.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
int SpectrumCmd::action(int argc, char **argv)
{
Wavefunction& wf = s->wf;
// Check that only the k=0 point is used
if ( wf.nkp()>1 || length(wf.kpoint(0)) != 0.0 )
{
if ( ui->onpe0() )
{
cout << " SpectrumCmd::action: spectrum only implemented at\n"
<< " the Gamma point (k=0)" << endl;
}
return 1;
}
if ( !( argc == 2 || argc == 3 || argc == 5 ) )
{
if ( ui->onpe0() )
{
cout << " use: spectrum [emin emax] [width] filename" << endl;
}
return 1;
}
// Compute eigenvalues using the current wave function wf
Wavefunction dwf(wf);
ChargeDensity cd(wf);
cd.update_density();
EnergyFunctional ef(*s,cd);
const bool compute_stress = false;
ef.update_vhxc(compute_stress);
const bool compute_forces = false;
const bool compute_hpsi = true;
valarray<double> sigma_eks;
vector<vector<double> > fion;
ef.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
const bool compute_eigvec = true;
wf.diag(dwf,compute_eigvec);
if ( ui->onpe0() )
{
cout << "<eigenset>" << endl;
// print eigenvalues
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const int nst = wf.sd(ispin,ikp)->nst();
const double eVolt = 2.0 * 13.6058;
cout << " <eigenvalues spin=\"" << ispin
<< "\" kpoint=\""
<< setprecision(8)
<< wf.sd(ispin,ikp)->kpoint()
<< "\" weight=\""
<< setprecision(8)
<< wf.weight(ikp)
<< "\" n=\"" << nst << "\">" << endl;
for ( int i = 0; i < nst; i++ )
{
cout << setw(12) << setprecision(5)
<< wf.sd(ispin,ikp)->eig(i)*eVolt;
if ( i%5 == 4 ) cout << endl;
}
if ( nst%5 != 0 ) cout << endl;
cout << " </eigenvalues>" << endl;
}
}
cout << "</eigenset>" << endl;
}
const double eVolt = 2.0 * 13.6058;
const UnitCell& cell = wf.cell();
// emin, emax: bounds of plot in eV
// de: energy spacing of plot values in eV (fixed at 0.01)
// width: gaussian width of convolution in eV (default 0.05)
const double de = 0.01;
double emin = 0.0, emax = 0.0, width = 0.05, erange = 0.0;
const char *spfilename;
// spectrum file
if ( argc == 2 )
{
spfilename = argv[1];
}
// spectrum width file
if ( argc == 3 )
{
width = atof(argv[1]);
spfilename = argv[2];
}
// spectrum emin emax de width file
if ( argc == 5 )
{
emin = atof(argv[1]);
emax = atof(argv[2]);
width = atof(argv[3]);
spfilename = argv[4];
erange = emax - emin + 3 * width;
if ( emax <= emin )
{
if ( ui->onpe0() )
{
cout << " SpectrumCmd::action: emax must be larger than emin" << endl;
}
return 1;
}
}
const int nspin = wf.nspin();
vector<vector<double> > sp(nspin);
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
SlaterDet& sd = *(wf.sd(ispin,0));
const int n = sd.nst();
if ( erange == 0.0 )
erange = eVolt * ( sd.eig(n-1) - sd.eig(0) ) + 3 * width;
const int np = erange / de;
sp[ispin].resize(np);
fill(sp[ispin].begin(),sp[ispin].end(),0.0);
MLWFTransform* mlwft = new MLWFTransform(sd);
mlwft->update();
if ( ui->onpe0() )
cout << "<dipole_matrix_elements ispin=\"" << ispin+1 << "\">" << endl;
const DoubleMatrix *amat[6];
for ( int i = 0; i < 6; i++ )
amat[i] = mlwft->a(i);
const double *c0 = amat[0]->cvalptr();
const double *s0 = amat[1]->cvalptr();
const double *c1 = amat[2]->cvalptr();
const double *s1 = amat[3]->cvalptr();
const double *c2 = amat[4]->cvalptr();
const double *s2 = amat[5]->cvalptr();
const Context& ctxt = amat[0]->context();
// loop over global indices (i,j)
for ( int i = 0; i < n; i++ )
{
for ( int j = i+1; j < n; j++ )
{
// if i occupied and j empty
if ( sd.occ(i) > 0.0 && sd.occ(j) == 0.0 )
{
const double delta_e = eVolt * ( sd.eig(j) - sd.eig(i) );
// dipole transition strength w
double w = 0.0;
// if element (i,j) is located on the current task,
// compute local indices (iloc,jloc)
const int pr = amat[0]->pr(i);
const int pc = amat[0]->pc(j);
if ( pr == ctxt.myrow() && pc == ctxt.mycol() )
{
const int iloc = amat[0]->l(i) * amat[0]->mb() + amat[0]->x(i);
const int jloc = amat[0]->m(j) * amat[0]->nb() + amat[0]->y(j);
const int mloc = amat[0]->mloc();
// position in local array is k = iloc+mloc*jloc
const int k = iloc + mloc * jloc;
const double fac = 0.5 * M_1_PI;
double c[3] = { fac*c0[k], fac*c1[k], fac*c2[k] };
double s[3] = { fac*s0[k], fac*s1[k], fac*s2[k] };
// cc, ss: matrix elements in cartesian coordinates
double cc[3], ss[3];
cell.vecmult3x3(cell.amat(),c,cc);
cell.vecmult3x3(cell.amat(),s,ss);
w = cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2]+
ss[0]*ss[0]+ss[1]*ss[1]+ss[2]*ss[2];
// add contribution to the absorption spectrum
for ( int ie = 0; ie < np; ie++ )
{
const double t = ( emin + ie * de - delta_e ) / width;
sp[ispin][ie] += w * width * sqrt(M_PI) * exp(-t*t);
}
// only send if not on pe 0
if ( !ui->onpe0() )
ctxt.dsend(1,1,&w,1,0,0);
} // if pr,pc
// receive if on pe 0 and element was sent from another pe
if ( ui->onpe0() && !( pr==0 && pc==0 ) )
ctxt.drecv(1,1,&w,1,pr,pc);
if ( ui->onpe0() )
{
cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::right, ios::adjustfield);
cout << " <dipole i=\"" << i+1 << "\" j=\"" << j+1 << "\"> ";
cout << setprecision(6) << delta_e << " "
<< w << " </dipole>" << endl;
}
} // if i occupied and j empty
} // for j
} // for i
// sum contributions to sp from all tasks
ctxt.dsum(np,1,&sp[ispin][0],np);
if ( ui->onpe0() )
cout << "</dipole_matrix_elements>" << endl;
delete mlwft;
} // ispin
// write spectrum to file spectrum.dat
if ( ui->onpe0() )
{
ofstream spfile(spfilename);
for ( int ispin = 0; ispin < nspin; ispin++ )
{
spfile << "# spectrum ispin=" << ispin+1
<< " width=" << width << endl;
for ( int ie = 0; ie < sp[ispin].size(); ie++ )
spfile << emin + ie * de << " " << sp[ispin][ie] << endl;
spfile << endl << endl;
}
spfile.close();
}
return 0;
}
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// SpectrumCmd.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef SPECTRUMCMD_H
#define SPECTRUMCMD_H
#include <iostream>
#include <stdlib.h>
#include <string>
#include "UserInterface.h"
class Sample;
class SpectrumCmd : public Cmd
{
private:
public:
Sample *s;
SpectrumCmd(Sample *sample) : s(sample) {};
const char *name(void) const { return "spectrum"; }
const char *help_msg(void) const
{
return
"\n spectrum\n\n"
" syntax: spectrum filename\n"
" spectrum width filename\n"
" spectrum emin emax width filename\n\n"
" The spectrum command computes the dipole transition strengths\n"
" between occupied and empty orbitals. It computes Kohn-Sham eigenvalues\n"
" and eigenfunctions of the current wave function using the the current\n"
" value of the xc variable. The corresponding absorption\n"
" spectrum is written on an output file after convolution with a\n"
" gaussian function.\n"
" emin, emax: energy range (optional)\n"
" width : width of the gaussian (optional) (default 0.05 eV)\n"
" filename : output file name\n"
" If emin and emax are not given, the energy range includes\n"
" all possible transitions between occupied and empty orbitals.\n"
" All energy parameters must be given in eV.\n\n";
}
int action(int argc, char **argv);
SpectrumCmd();
};
#endif
--------------------------------------------------------------------------------
rel1_69_0
--------------------------------------------------------------------------------
beb8bee Cleanup trailing spaces in util scripts
3137cbd Add scripts qbox_reduce.sh, qbox_move_subsample.py
8d434af Fix spectrum file for nspin==2
7a00811 Update spectrum command
6ba92e5 Update spectrum command help msg
80c343d Initialize eig with zero in SlaterDet
4031b74 Modify spectrum command parameters
539499d write spectrum on file, fix parallel SpectrumCmd
54058fe Add spectrum cmd to Makefile and qb.C
fc34625 Add spectrum command
--------------------------------------------------------------------------------
rel1_68_4
--------------------------------------------------------------------------------
d9bac45 Use syevd/heevd in JD stepper, workaround lwork size calculation
......
......@@ -70,6 +70,7 @@ using namespace std;
#include "SetCmd.h"
#include "SetVelocityCmd.h"
#include "SpeciesCmd.h"
#include "SpectrumCmd.h"
#include "StatusCmd.h"
#include "StrainCmd.h"
#include "TorsionCmd.h"
......@@ -264,6 +265,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new SetCmd(s));
ui.addCmd(new SetVelocityCmd(s));
ui.addCmd(new SpeciesCmd(s));
ui.addCmd(new SpectrumCmd(s));
ui.addCmd(new StatusCmd(s));
ui.addCmd(new StrainCmd(s));
ui.addCmd(new TorsionCmd(s));
......
......@@ -19,5 +19,5 @@
#include "release.h"
std::string release(void)
{
return std::string("1.68.4");
return std::string("rel1_69_0");
}
......@@ -4,5 +4,5 @@
# use: get_atomset sample.xml
#
nlines=$(grep /atomset -m 1 -n $1 | cut -f1 -d: - )
head -$nlines $1
head -$nlines $1
echo "</fpmd:sample>"
......@@ -36,7 +36,7 @@ de = (emax - emin)/(ndos-1)
# f(x) = 1/(sqrt(pi)*width) * exp(-(x/width)^2 )
def gauss(x, width):
return (1.0/(math.sqrt(math.pi)*width)) * math.exp(-(x/width)**2)
# Qbox output handler to extract and process data
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
......
#!/usr/bin/python
# Copyright 2018 The Regents of the University of California
# This file is part of Qbox
#
# qbox_move_subsample.py: create a qbox input file moving atoms to the positions
# of iterations of a simulation subsampled every "interval" number of steps
# Each set of moves is followed by a command cmd given as an argument
#
# use: qbox_move_subsample.py interval cmd {file|URL}
import os.path
import xml.sax
import sys
import urllib2
def usage():
print "use: ",sys.argv[0]," interval cmd {file|URL}"
sys.exit()
argc=len(sys.argv)
if ( argc != 4 ):
usage()
interval = int(sys.argv[1])
cmd = sys.argv[2]
input_source = sys.argv[3]
# Qbox output handler to extract and process data
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
self.step = 0
self.inAtomset = 0
self.inAtom = 0
self.inPosition = 0
def startElement(self, name, attributes):
if name == "atomset":
self.tau=[]
self.atomname=[]
self.inAtomset = 1
elif (name == "unit_cell") & self.inAtomset:
self.cell_a = attributes["a"]
self.cell_b = attributes["b"]
self.cell_c = attributes["c"]
elif (name == "atom") & self.inAtomset:
self.atomname.append(attributes["name"])
self.inAtom = 1
elif (name == "position") & self.inAtom:
self.buffer = ""
self.inPosition = 1
def characters(self, data):
if self.inPosition:
self.buffer += data
def endElement(self, name):
if (name == "atom") and self.inAtomset:
self.inAtom = 0
if (name == "position") & self.inAtom:
pos = self.buffer.split()
self.tau.append([pos[0],pos[1],pos[2]])
self.inPosition = 0
elif name == "atomset":
self.step += 1
if ( self.step % interval == 0 ):
print "#",input_source,"iteration",self.step
avec = self.cell_a.split()
bvec = self.cell_b.split()
cvec = self.cell_c.split()
print "set cell ",avec[0],avec[1],avec[2],\
bvec[0],bvec[1],bvec[2],\
cvec[0],cvec[1],cvec[2]
for i in range(len(self.tau)):
print "move ",self.atomname[i]," to ",\
self.tau[i][0],self.tau[i][1],self.tau[i][2]
print cmd
self.inAtomset = 0
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
# test if input_source is a local file
# if not, process as a URL
if ( os.path.isfile(input_source) ):
file = open(input_source)
s = file.read(8192)
while ( s !="" ):
parser.feed(s)
s = file.read(8192)
file.close()
else:
# attempt to open as a URL
try:
f = urllib2.urlopen(input_source)
s = f.read(8192)
while ( s !="" ):
parser.feed(s)
s = f.read(8192)
f.close()
except (ValueError,urllib2.HTTPError) as e:
print e
sys.exit()
parser.reset()
#!/bin/bash
# reduce qbox restart file to atomset file
# use: qbox_reduce.sh file.xml [file.xml ..]
for f in ${*}
do
name=${f%.xml}
atomset_name=${name}_atomset
echo $name.xml "->" $atomset_name.xml
get_atomset $f > $atomset_name.xml
rm $f
done
......@@ -73,9 +73,9 @@ for h in hlist:
print "# current ",h[1]," at ", h[3],h[4],h[5]
print "# nearest O is at ", nearest_o[3],nearest_o[4],nearest_o[5]
print "# move ",h[1]," by ", sx_min, sy_min, sz_min
h[3] += sx_min
h[4] += sy_min
h[5] += sz_min
h[3] += sx_min
h[4] += sy_min
h[5] += sz_min
for o in olist:
print o[0],o[1],o[2],'%10.5f'%o[3],'%10.5f'%o[4],'%10.5f'%o[5]
......
......@@ -7,7 +7,7 @@
# use: qbox_species_temp.sh species_name file [file ...]
#
if (($#<2))
then echo " use: qbox_species_temp.sh species file [file ...]"
then echo " use: qbox_species_temp.sh species file [file ...]"
exit 1
fi
species=$1
......
#!/bin/bash
# qbox_translate: translate all atoms
# qbox_translate: translate all atoms
# use: qbox_replicate cell.sys dx dy dz > newcell.sys
#
if (( $# != 4 ))
......
......@@ -3,7 +3,7 @@
# This file is part of Qbox
#
# qbox_xyz.py: extract sets of atomic positions in xyz format
# from a Qbox output file or from a Qbox sample file using SAX
# from a Qbox output file or from a Qbox sample file using SAX
# incremental parsing
#
# use: qbox_xyz.py [-first] {file|URL}
......
......@@ -45,10 +45,10 @@ class AtomSet:
class Sample:
def __init__(self):
self.atoms = AtomSet()
# The following handler processes the <atomset> element of
# an XML document and updates the AtomSet data of the Sample
# If multiple instances of <atomset> are found, the
# If multiple instances of <atomset> are found, the
# handler overwrites the AtomSet data
class QSOAtomSetHandler(xml.sax.handler.ContentHandler):
def __init__(self,sample):
......@@ -73,7 +73,7 @@ class QSOAtomSetHandler(xml.sax.handler.ContentHandler):
self.s.atoms.cell.a = attributes["a"]
self.s.atoms.cell.b = attributes["b"]
self.s.atoms.cell.c = attributes["c"]
elif (name == "species"):
elif (name == "species"):
self.inSpecies = True
self.species_name = "species_name"
if "name" in attributes:
......@@ -83,9 +83,9 @@ class QSOAtomSetHandler(xml.sax.handler.ContentHandler):
self.species_href = attributes["href"]
self.species_symbol = self.species_name+"_symbol"
self.species_atomic_number = self.species_name+"_atomic_number"
self.species_mass = self.species_name+"_mass"
self.species_mass = self.species_name+"_mass"
elif (name == "atom") and self.inAtomSet:
self.inAtom = True
self.inAtom = True
self.atom_name = attributes["name"]
self.atom_species = attributes["species"]
sp = Species(self.species_name,self.species_href,self.species_symbol,
......@@ -124,7 +124,7 @@ class QSOAtomSetHandler(xml.sax.handler.ContentHandler):
x = float(pos[0])
y = float(pos[1])
z = float(pos[2])
self.atom_position = [x,y,z]
self.atom_position = [x,y,z]
self.inPosition = False
if (name == "velocity") and self.inAtom:
vel = self.buffer.split()
......
#!/usr/bin/python
# Convert <atomset> elements from quantum-simulation.org (QSO) format
# Convert <atomset> elements from quantum-simulation.org (QSO) format
# to Qbox input file
# use: qso2qbox.py [-last] {file|URL}
# Default: only the first <atomset> element is processed
......@@ -10,7 +10,7 @@ import os.path
import xml.sax
import sys
import urllib2
import datetime
import datetime
def usage():
print "use: ",sys.argv[0]," [-last] {file|URL}"
......
......@@ -18,7 +18,7 @@ xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0">
</xsl:template>
<xsl:template match="atom">
<xsl:text>move </xsl:text>
<xsl:value-of select="@name"/>
<xsl:value-of select="@name"/>
<xsl:text> to </xsl:text>
<xsl:value-of select="position"/> <xsl:text>
</xsl:text>
......
......@@ -32,8 +32,8 @@ xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0">
<xsl:text>atom </xsl:text>
<xsl:value-of select="@name"/> <xsl:text> </xsl:text>
<xsl:value-of select="@species"/> <xsl:text> </xsl:text>
<xsl:value-of select="position"/> <xsl:text> </xsl:text>
<xsl:value-of select="velocity"/> <xsl:text>