Commit 5d871973 by Francois Gygi

upf2qso translation program


git-svn-id: http://qboxcode.org/svn/qb/trunk@1434 cba15fb0-1239-40c8-b417-11db7ca47a34
parent a48d8450
upf2qso.C is a program that converts norm-conserving pseudopotentials
written in the UPF format (see http://www.pwscf.org) to the
www.quantum-simulation.org (QSO) format.
upf2qso only works for norm-conserving potentials. Ultrasoft potentials
are not supported
Test:
$ src/upf2qso 5.0 < data/C.pz-vbc.UPF > C.pz-vbc.xml
$ src/upf2qso 8.0 < data/Si.pz-vbc.UPF > Si.pz-vbc.xml
$ src/upf2qso 8.0 < data/30-Zn.LDA.fhi.UPF > Zn_FHI_LDA.xml
Additional files v.dat, vlin.dat, upf.dat contain various functions in
gnuplot compatible format
This program was written by F. Gygi
Debugging of version 1.2 was provided by Stefan Wippermann
Version 1.4, 2014-02-24. Fixed bug in 1.3 in calc of lmax
CXXFLAGS=-g
upf2qso: upf2qso.o PeriodicTable.o spline.o
$(CXX) -o $@ $^
clean:
rm -f *.o upf2qso
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PeriodicTable.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PeriodicTable.C,v 1.3 2009/02/14 22:47:30 fgygi Exp $
#include "PeriodicTable.h"
#include <cassert>
////////////////////////////////////////////////////////////////////////////////
int PeriodicTable::z(string symbol) const
{
map<string,int>::const_iterator i = zmap.find(symbol);
assert( i != zmap.end() );
return (*i).second;
}
////////////////////////////////////////////////////////////////////////////////
string PeriodicTable::symbol(int z) const
{
assert(z>0 && z<=ptable.size());
return ptable[z-1].symbol;
}
////////////////////////////////////////////////////////////////////////////////
string PeriodicTable::configuration(int z) const
{
assert(z>0 && z<=ptable.size());
return ptable[z-1].config;
}
////////////////////////////////////////////////////////////////////////////////
string PeriodicTable::configuration(string symbol) const
{
return ptable[z(symbol)-1].config;
}
////////////////////////////////////////////////////////////////////////////////
double PeriodicTable::mass(int z) const
{
assert(z>0 && z<=ptable.size());
return ptable[z-1].mass;
}
////////////////////////////////////////////////////////////////////////////////
double PeriodicTable::mass(string symbol) const
{
return ptable[z(symbol)-1].mass;
}
////////////////////////////////////////////////////////////////////////////////
int PeriodicTable::size(void) const
{
return ptable.size();
}
////////////////////////////////////////////////////////////////////////////////
PeriodicTable::PeriodicTable(void)
{
ptable.push_back(Element(1,"H","1s1",1.00794));
ptable.push_back(Element(2,"He","1s2",4.00260));
ptable.push_back(Element(3, "Li","1s2 2s1", 6.941));
ptable.push_back(Element(4, "Be","1s2 2s2", 9.01218));
ptable.push_back(Element(5, "B", "1s2 2s2 2p1",10.811));
ptable.push_back(Element(6, "C", "1s2 2s2 2p2",12.0107));
ptable.push_back(Element(7, "N", "1s2 2s2 2p3",14.00674));
ptable.push_back(Element(8, "O", "1s2 2s2 2p4",15.9994));
ptable.push_back(Element(9, "F", "1s2 2s2 2p5",18.9884));
ptable.push_back(Element(10,"Ne","1s2 2s2 2p6",20.1797));
ptable.push_back(Element(11,"Na","[Ne] 3s1", 22.98977));
ptable.push_back(Element(12,"Mg","[Ne] 3s2", 24.3050));
ptable.push_back(Element(13,"Al","[Ne] 3s2 3p1",26.98154));
ptable.push_back(Element(14,"Si","[Ne] 3s2 3p2",28.0855));
ptable.push_back(Element(15,"P", "[Ne] 3s2 3p3",30.97376));
ptable.push_back(Element(16,"S", "[Ne] 3s2 3p4",32.066));
ptable.push_back(Element(17,"Cl","[Ne] 3s2 3p5",35.4527));
ptable.push_back(Element(18,"Ar","[Ne] 3s2 3p6",39.948));
ptable.push_back(Element(19,"K", "[Ar] 4s1",39.0983));
ptable.push_back(Element(20,"Ca","[Ar] 4s2",40.078));
ptable.push_back(Element(21,"Sc","[Ar] 3d1 4s2",44.95591));
ptable.push_back(Element(22,"Ti","[Ar] 3d2 4s2",47.867));
ptable.push_back(Element(23,"V", "[Ar] 3d3 4s2",50.9415));
ptable.push_back(Element(24,"Cr","[Ar] 3d5 4s1",51.9961));
ptable.push_back(Element(25,"Mn","[Ar] 3d5 4s2",54.93805));
ptable.push_back(Element(26,"Fe","[Ar] 3d6 4s2",55.845));
ptable.push_back(Element(27,"Co","[Ar] 3d7 4s2",58.9332));
ptable.push_back(Element(28,"Ni","[Ar] 3d8 4s2",58.6934));
ptable.push_back(Element(29,"Cu","[Ar] 3d10 4s1",63.546));
ptable.push_back(Element(30,"Zn","[Ar] 3d10 4s2",65.39));
ptable.push_back(Element(31,"Ga","[Ar] 3d10 4s2 4p1",69.723));
ptable.push_back(Element(32,"Ge","[Ar] 3d10 4s2 4p2",72.61));
ptable.push_back(Element(33,"As","[Ar] 3d10 4s2 4p3",74.9216));
ptable.push_back(Element(34,"Se","[Ar] 3d10 4s2 4p4",78.96));
ptable.push_back(Element(35,"Br","[Ar] 3d10 4s2 4p5",79.904));
ptable.push_back(Element(36,"Kr","[Ar] 3d10 4s2 4p6",83.80));
ptable.push_back(Element(37,"Rb","[Kr] 5s1",85.4678));
ptable.push_back(Element(38,"Sr","[Kr] 5s2",87.62));
ptable.push_back(Element(39,"Y" ,"[Kr] 4d1 5s2",88.90585));
ptable.push_back(Element(40,"Zr","[Kr] 4d2 5s2",91.224));
ptable.push_back(Element(41,"Nb","[Kr] 4d4 5s1",92.90638));
ptable.push_back(Element(42,"Mo","[Kr] 4d5 5s1",95.94));
ptable.push_back(Element(43,"Tc","[Kr] 4d5 5s2",98.0));
ptable.push_back(Element(44,"Ru","[Kr] 4d7 5s1",101.07));
ptable.push_back(Element(45,"Rh","[Kr] 4d8 5s1",102.9055));
ptable.push_back(Element(46,"Pd","[Kr] 4d10",106.42));
ptable.push_back(Element(47,"Ag","[Kr] 4d10 5s1",107.8682));
ptable.push_back(Element(48,"Cd","[Kr] 4d10 5s2",112.411));
ptable.push_back(Element(49,"In","[Kr] 4d10 5s2 5p1",114.818));
ptable.push_back(Element(50,"Sn","[Kr] 4d10 5s2 5p2",118.710));
ptable.push_back(Element(51,"Sb","[Kr] 4d10 5s2 5p3",121.760));
ptable.push_back(Element(52,"Te","[Kr] 4d10 5s2 5p4",127.60));
ptable.push_back(Element(53,"I" ,"[Kr] 4d10 5s2 5p5",126.90447));
ptable.push_back(Element(54,"Xe","[Kr] 4d10 5s2 5p6",131.29));
ptable.push_back(Element(55,"Cs","[Xe] 6s1",132.90545));
ptable.push_back(Element(56,"Ba","[Xe] 6s2",137.327));
ptable.push_back(Element(57,"La","[Xe] 5d1 6s2",138.9055));
ptable.push_back(Element(58,"Ce","[Xe] 4f1 5d1 6s2",140.116));
ptable.push_back(Element(59,"Pr","[Xe] 4f3 6s2",140.90765));
ptable.push_back(Element(60,"Nd","[Xe] 4f4 6s2",144.24));
ptable.push_back(Element(61,"Pm","[Xe] 4f5 6s2",145.0));
ptable.push_back(Element(62,"Sm","[Xe] 4f6 6s2",150.36));
ptable.push_back(Element(63,"Eu","[Xe] 4f7 6s2",151.964));
ptable.push_back(Element(64,"Gd","[Xe] 4f7 5d1 6s2",157.25));
ptable.push_back(Element(65,"Tb","[Xe] 4f9 6s2",158.92534));
ptable.push_back(Element(66,"Dy","[Xe] 4f10 6s2",162.50));
ptable.push_back(Element(67,"Ho","[Xe] 4f11 6s2",164.93032));
ptable.push_back(Element(68,"Er","[Xe] 4f12 6s2",167.26));
ptable.push_back(Element(69,"Tm","[Xe] 4f13 6s2",168.93421));
ptable.push_back(Element(70,"Yb","[Xe] 4f14 6s2",173.04));
ptable.push_back(Element(71,"Lu","[Xe] 4f14 5d1 6s2",174.967));
ptable.push_back(Element(72,"Hf","[Xe] 4f14 5d2 6s2",178.49));
ptable.push_back(Element(73,"Ta","[Xe] 4f14 5d3 6s2",180.9479));
ptable.push_back(Element(74,"W" ,"[Xe] 4f14 5d4 6s2",183.84));
ptable.push_back(Element(75,"Re","[Xe] 4f14 5d5 6s2",186.207));
ptable.push_back(Element(76,"Os","[Xe] 4f14 5d6 6s2",190.23));
ptable.push_back(Element(77,"Ir","[Xe] 4f14 5d7 6s2",192.217));
ptable.push_back(Element(78,"Pt","[Xe] 4f14 5d9 6s1",195.078));
ptable.push_back(Element(79,"Au","[Xe] 4f14 5d10 6s1",196.96655));
ptable.push_back(Element(80,"Hg","[Xe] 4f14 5d10 6s2",200.59));
ptable.push_back(Element(81,"Tl","[Xe] 4f14 5d10 6s2 6p1",204.3833));
ptable.push_back(Element(82,"Pb","[Xe] 4f14 5d10 6s2 6p2",207.2));
ptable.push_back(Element(83,"Bi","[Xe] 4f14 5d10 6s2 6p3",208.98038));
ptable.push_back(Element(84,"Po","[Xe] 4f14 5d10 6s2 6p4",209.0));
ptable.push_back(Element(85,"At","[Xe] 4f14 5d10 6s2 6p5",210.0));
ptable.push_back(Element(86,"Rn","[Xe] 4f14 5d10 6s2 6p6",222.0));
ptable.push_back(Element(87,"Fr","[Rn] 7s1",223.0));
ptable.push_back(Element(88,"Ra","[Rn] 7s2",226.0));
ptable.push_back(Element(89,"Ac","[Rn] 6d1 7s2",227.0));
ptable.push_back(Element(90,"Th","[Rn] 6d2 7s2",232.0381));
ptable.push_back(Element(91,"Pa","[Rn] 5f2 6d1 7s2",231.03588));
ptable.push_back(Element(92,"U" ,"[Rn] 5f3 6d1 7s2",238.0289));
ptable.push_back(Element(93,"Np","[Rn] 5f4 6d1 7s2",237.0));
ptable.push_back(Element(94,"Pu","[Rn] 5f5 6d1 7s2",244.0));
for ( int i = 0; i < ptable.size(); i++ )
zmap[ptable[i].symbol] = i+1;
}
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2009 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// PeriodicTable.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PeriodicTable.h,v 1.2 2009/02/14 22:47:30 fgygi Exp $
#ifndef PERIODICTABLE_H
#define PERIODICTABLE_H
#include <map>
#include <vector>
#include <string>
using namespace std;
struct Element
{
int z;
string symbol;
string config;
double mass;
Element(int zz, string s, string c, double m) : z(zz), symbol(s), config(c),
mass(m) {}
};
class PeriodicTable
{
private:
vector<Element> ptable;
map<string,int> zmap;
public:
PeriodicTable(void);
int z(string symbol) const;
string symbol(int zval) const;
string configuration(int zval) const;
string configuration(string symbol) const;
double mass(int zval) const;
double mass(string symbol) const;
int size(void) const;
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
///////////////////////////////////////////////////////////////////////////////
//
// spline.C
//
///////////////////////////////////////////////////////////////////////////////
// $Id: spline.C,v 1.5 2008/09/08 15:56:20 fgygi Exp $
#include "spline.h"
#include <cassert>
void tridsolve(int n, double* d, double* e, double* f, double* x)
{
// solve the tridiagonal system Ax=b
// d[i] = a(i,i)
// e[i] = a(i,i+1) (superdiagonal of A, e[n-1] not defined)
// f[i] = a(i,i-1) (subdiagonal of A, f[0] not defined)
// x[i] = right-hand side b as input
// x[i] = solution as output
for ( int i = 1; i < n; i++ )
{
f[i] /= d[i-1];
d[i] -= f[i]*e[i-1];
}
for ( int i = 1; i < n; i++ )
x[i] -= f[i]*x[i-1];
x[n-1] /= d[n-1];
for ( int i = n-2; i >= 0; i-- )
x[i] = (x[i]-e[i]*x[i+1])/d[i];
}
void spline(int n, double *x, double *y, double yp_left, double yp_right,
int bcnat_left, int bcnat_right, double *y2)
{
const double third = 1.0/3.0;
const double sixth = 1.0/6.0;
double *d = new double[n];
double *e = new double[n];
double *f = new double[n];
if ( bcnat_left == 0 )
{
// use derivative yp_left at x[0]
const double h = x[1]-x[0];
assert(h>0.0);
d[0] = third*h;
e[0] = sixth*h;
f[0] = 0.0;
y2[0] = (y[1]-y[0])/h - yp_left;
}
else
{
// use natural spline at x[0]
d[0] = 1.0;
e[0] = 0.0;
f[0] = 0.0;
y2[0] = 0.0;
}
if ( bcnat_right == 0 )
{
// use derivative yp_right at x[n-1]
const double h = x[n-1]-x[n-2];
assert(h>0.0);
d[n-1] = third*h;
e[n-1] = 0.0;
f[n-1] = sixth*h;
y2[n-1] = yp_right - (y[n-1]-y[n-2])/h;
}
else
{
// use natural spline at x[n-1]
d[n-1] = 1.0;
e[n-1] = 0.0;
f[n-1] = 0.0;
y2[n-1] = 0.0;
}
// tridiagonal matrix
for ( int i = 1; i < n-1; i++ )
{
const double hp = x[i+1]-x[i];
const double hm = x[i]-x[i-1];
assert(hp>0.0);
assert(hm>0.0);
d[i] = third * (hp+hm);
e[i] = sixth * hp;
f[i] = sixth * hm;
y2[i] = (y[i+1]-y[i])/hp - (y[i]-y[i-1])/hm;
}
tridsolve(n,d,e,f,y2);
delete [] d;
delete [] e;
delete [] f;
}
void splint (int n, double *xa, double *ya, double *y2a, double x, double *y)
{
int k;
double a,b,h;
int kl = 0;
int kh = n-1;
while ( kh - kl > 1 )
{
k = ( kh + kl ) / 2;
if ( xa[k] > x )
kh = k;
else
kl = k;
}
h = xa[kh] - xa[kl];
assert ( h > 0.0 );
a = ( xa[kh] - x ) / h;
b = ( x - xa[kl] ) / h;
*y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
}
void splintd (int n, double *xa, double *ya, double *y2a,
double x, double *y, double *dy)
{
int k;
double a,b,h;
int kl = 0;
int kh = n-1;
while ( kh - kl > 1 )
{
k = ( kh + kl ) / 2;
if ( xa[k] > x )
kh = k;
else
kl = k;
}
h = xa[kh] - xa[kl];
assert ( h > 0.0 );
a = ( xa[kh] - x ) / h;
b = ( x - xa[kl] ) / h;
*y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
*dy = ( ya[kh] - ya[kl] ) / h +
h * ( ( (1.0/6.0) - 0.5 * a * a ) * y2a[kl] +
( 0.5 * b * b - (1.0/6.0) ) * y2a[kh] );
}
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
///////////////////////////////////////////////////////////////////////////////
//
// spline.h
//
///////////////////////////////////////////////////////////////////////////////
// $Id: spline.h,v 1.5 2008/09/08 15:56:20 fgygi Exp $
void spline(int n, double *x, double *y, double yp_left, double yp_right,
int bcnat_left, int bcnat_right, double *y2);
void splint (int n, double *xa, double *ya, double *y2a, double x, double *y);
void splintd (int n, double *xa, double *ya, double *y2a,
double x, double *y, double *dy);
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