Commit f6060b32 by Francois Gygi

Modified ResponseCmd to use xml vext file

ExternalPotential uses a Function3d object to read a vext XML file.
ResponseCmd option modified: -io not required, default xml.
-nx,-ny,-nz not used, size read from either cube or xml file.
parent d9bc9a40
......@@ -29,6 +29,7 @@ using namespace std;
#include "Basis.h"
#include "ExternalPotential.h"
#include "FourierTransform.h"
#include "Function3d.h"
#include "Base64Transcoder.h"
bool abs_compare(const double &a, const double &b)
......@@ -49,10 +50,8 @@ void ExternalPotential::update(const ChargeDensity& cd)
Timer tm_read_vext, tm_comm_vext;
double time, tmin, tmax;
// In cube mode and base64_serial mode, the whole external potential is
// In cube mode and xml mode, the whole external potential is
// read by all processors on row 1 and stored in vext_read
// In base64_parallel mode, each processor on column 1 read part of external
// potential and store in vext_read_loc
vector<double> vext_read, vext_read_loc;
tm_read_vext.start();
......@@ -61,7 +60,7 @@ void ExternalPotential::update(const ChargeDensity& cd)
// read cube file, n_'s are determined by cube file
if ( myrow == 0 )
{
ifstream vfile(filename_.c_str(), std::ios::in);
ifstream vfile(filename_.c_str());
if (!vfile)
{
if (mycol == 0)
......@@ -80,10 +79,9 @@ void ExternalPotential::update(const ChargeDensity& cd)
vfile >> n_[i];
getline(vfile, tmpstr);
}
n012_ = n_[0] * n_[1] * n_[2];
for (int i = 0; i < natom; i++)
getline(vfile, tmpstr); // skip atom coordinates
vext_read.resize(n012_);
vext_read.resize(n_[0] * n_[1] * n_[2]);
for (int nx = 0; nx < n_[0]; nx++)
for (int ny = 0; ny < n_[1]; ny++)
for (int nz = 0; nz < n_[2]; nz++)
......@@ -94,35 +92,19 @@ void ExternalPotential::update(const ChargeDensity& cd)
vfile.close();
}
MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
MPI_Bcast(&n012_,1,MPI_INT,0,vcomm);
}
else if (io_ == "base64_serial")
else if (io_ == "xml")
{
if (myrow == 0)
{
ifstream vfile(filename_.c_str(), std::ios::in | std::ios::binary);
if (!vfile)
{
if (mycol == 0)
cout << " ExternalPotential::update: file not found: "
<< filename_ << endl;
ctxt->abort(1);
}
Base64Transcoder xcdr;
int nchars = xcdr.nchars(n012_ * sizeof(double));
char *rbuf = new char[nchars];
vfile.seekg(0, ios::end);
assert(vfile.tellg() == nchars);
vfile.seekg(0, ios::beg);
vfile.read(rbuf, nchars);
// vector<char> rbuf(istreambuf_iterator<char>(vfile),
// istreambuf_iterator<char>());
// assert(rbuf.size() == nchars);
vext_read.resize(n012_);
xcdr.decode(nchars, rbuf, (byte *) (&vext_read[0]));
vfile.close();
delete [] rbuf;
Function3d f;
f.read(filename_);
vext_read = f.val;
n_[0] = f.nx;
n_[1] = f.ny;
n_[2] = f.nz;
}
MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
}
tm_read_vext.stop();
// at this point, all processes have correct n_ regardless of io mode
......@@ -160,105 +142,21 @@ void ExternalPotential::update(const ChargeDensity& cd)
FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
vext_r_.resize(ft2.np012loc());
if ( io_== "base64_parallel" )
{
// base64_parallel mode: all processors on column 0 read
// in parallel, then bcast to other columns
int np012 = ft1.np012();
int np012loc = ft1.np012loc();
int lastproc = nprow - 1;
while ( lastproc >=0 && ft1.np2_loc(lastproc) == 0 ) lastproc --;
if ( mycol == 0 )
{
tm_read_vext.start();
Base64Transcoder xcdr;
int np012end;
MPI_Scan(&np012loc, &np012end, 1, MPI_INT, MPI_SUM, vcomm);
int np012start = np012end - np012loc;
np012end -= 1; // because index start from 0
// for simplicity, let the size of file (in bytes) read by each processor
// to be N*sizeof(double)*4, which correspond to N*3 float numbers
// denote each 3 floats as a group
int groupstart = np012start / 3;
int groupend = np012end / 3;
int offset = 4*sizeof(double)*groupstart;
int nfloats, nbytes, nchars;
if ( myrow < lastproc )
{
nfloats = 3 * (groupend - groupstart + 1);
}
else if ( myrow == lastproc )
{
nfloats = 3 * (groupend - groupstart) + (np012 - 1) % 3 + 1;
}
else
{
nfloats = 0;
}
if ( myrow < lastproc ) assert(nfloats >= vext_read_loc.size());
nbytes = nfloats * sizeof(double);
nchars = xcdr.nchars(nbytes);
char* rbuf = new char[nchars];
MPI_File fh;
MPI_Info info;
MPI_Info_create(&info);
int err;
err = MPI_File_open(vcomm, (char*) filename_.c_str(),
MPI_MODE_RDONLY , info, &fh);
if (err != 0)
cout << myrow << ": error in MPI_File_open: " << err << endl;
MPI_Status status;
err = MPI_File_read_at_all(fh,offset,(void*) &rbuf[0],nchars,
MPI_CHAR,&status);
if ( err != 0 )
cout << myrow << ": error in MPI_File_read_at_all: err=" << err << endl;
err = MPI_File_close(&fh);
if ( err != 0 )
cout << myrow << ": error in MPI_File_close: " << err << endl;
vector<double> tmpr(nfloats);
xcdr.decode(nchars, rbuf, (byte*) &tmpr[0]);
int istart = np012start % 3;
int iend = istart + np012loc - 1;
int j = 0;
for ( int i = istart; i <= iend; i++ )
vext_read_loc[j++] = tmpr[i];
delete [] rbuf;
tm_read_vext.stop();
} // if mycol == 0
tm_comm_vext.start();
// bcast vext_read_loc to other columns
if ( mycol == 0 )
ctxt->dbcast_send('r',np012loc,1,&vext_read_loc[0],np012loc);
else
ctxt->dbcast_recv('r',np012loc,1,&vext_read_loc[0],np012loc, myrow, 0);
tm_comm_vext.stop();
}
else
// xml mode or cube mode: processors on row 0 scatter
// vext to other rows
tm_comm_vext.start();
vector<int> scounts(nprow,0);
vector<int> sdispls(nprow,0);
int displ = 0;
for ( int iproc = 0; iproc < nprow; iproc++ )
{
// base64_serial mode or cube mode: processors on row 0 scatter
// vext to other rows
tm_comm_vext.start();
vector<int> scounts(nprow,0);
vector<int> sdispls(nprow,0);
int displ = 0;
for ( int iproc = 0; iproc < nprow; iproc++ )
{
sdispls[iproc] = displ;
scounts[iproc] = ft1.np012loc(iproc);
displ += ft1.np012loc(iproc);
}
MPI_Scatterv(&vext_read[0],&scounts[0],&sdispls[0],MPI_DOUBLE,
&vext_read_loc[0],ft1.np012loc(),MPI_DOUBLE,0,vcomm);
tm_comm_vext.stop();
sdispls[iproc] = displ;
scounts[iproc] = ft1.np012loc(iproc);
displ += ft1.np012loc(iproc);
}
MPI_Scatterv(&vext_read[0],&scounts[0],&sdispls[0],MPI_DOUBLE,
&vext_read_loc[0],ft1.np012loc(),MPI_DOUBLE,0,vcomm);
tm_comm_vext.stop();
// now vext_read_loc on all processors contains the correct portion of vext
// Fourier forward transform vext_read_loc to vext_g
......
......@@ -32,8 +32,7 @@ class ExternalPotential
Sample& s_;
int n_[3]; // real space grid size in 3 dimensions
// read from cube file in cube file mode,
// otherwise must be given in constructer
int n012_;
// otherwise must be given in constructor
double ecut_;
double magnitude_; // the magnitude of external potential, defined as
// the average of its largest 0.1% (absolute) values
......@@ -44,19 +43,11 @@ class ExternalPotential
public:
ExternalPotential(Sample& s, std::string name, std::string io="cube",
int nx=0, int ny=0, int nz=0):
ExternalPotential(Sample& s, std::string name, std::string io="xml"):
s_(s), filename_(name), ecut_(0.0), amplitude_(1.0), magnitude_(0.0),
io_(io)
{
assert( io_ == "cube" || io == "base64_serial" || io == "base64_parallel" );
if (io != "cube")
{
n_[0] = nx;
n_[1] = ny;
n_[2] = nz;
n012_ = n_[0] * n_[1] * n_[2];
}
assert( io_ == "cube" || io == "xml" );
}
~ExternalPotential() {}
......
......@@ -52,6 +52,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
uuid_str.o sampling.o CGOptimizer.o LineMinimizer.o \
ElectricEnthalpy.o PartialChargeCmd.o \
ExternalPotential.o ResponseCmd.o \
Function3d.o Function3dHandler.o \
$(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
$(EXEC): $(OBJECTS)
......
......@@ -37,9 +37,10 @@ using namespace std;
int ResponseCmd::action(int argc, char **argv)
{
// " syntax: response amplitude nitscf [nite]\n\n"
// " syntax: response -vext vext_file [-RPA|-IPA] [-amplitude a]
// " syntax: response -vext vext_file -drho drho_file
// [-RPA|-IPA] [-amplitude a]
// [-io iomode -nx nx -ny ny -nz nz]
// [-q qx qy qz] nitscf [nite]\n\n"
// nitscf [nite]\n\n"
if ( s->wf.nst() == 0 )
{
......@@ -72,9 +73,7 @@ int ResponseCmd::action(int argc, char **argv)
bool rpa = false;
bool ipa = false;
double amplitude = 0.0;
string io = "cube";
int nx, ny, nz;
nx = ny = nz = 0;
string io = "xml";
iarg++;
filename = argv[iarg];
......@@ -107,23 +106,23 @@ int ResponseCmd::action(int argc, char **argv)
{
iarg++;
io = argv[iarg];
assert( io == "cube" || io == "base64_serial" || io == "base64_parallel");
assert( io == "cube" || io == "xml" );
iarg++;
if ( io == "base64_serial" || io == "base64_parallel" )
{
assert(!strcmp(argv[iarg],"-nx"));
iarg++;
nx = atoi(argv[iarg]);
iarg++;
assert(!strcmp(argv[iarg],"-ny"));
iarg++;
ny = atoi(argv[iarg]);
iarg++;
assert(!strcmp(argv[iarg],"-nz"));
iarg++;
nz = atoi(argv[iarg]);
iarg++;
}
// if ( io == "xml" )
// {
// assert(!strcmp(argv[iarg],"-nx"));
// iarg++;
// nx = atoi(argv[iarg]);
// iarg++;
// assert(!strcmp(argv[iarg],"-ny"));
// iarg++;
// ny = atoi(argv[iarg]);
// iarg++;
// assert(!strcmp(argv[iarg],"-nz"));
// iarg++;
// nz = atoi(argv[iarg]);
// iarg++;
// }
}
nitscf = atoi(argv[iarg]);
......@@ -139,7 +138,7 @@ int ResponseCmd::action(int argc, char **argv)
return 1;
}
s->vext = new ExternalPotential(*s, filename, io, nx, ny, nz);
s->vext = new ExternalPotential(*s, filename, io);
s->vext->set_amplitude(amplitude);
responseVext(rpa, ipa, nitscf, nite, io);
delete s->vext;
......@@ -148,6 +147,7 @@ int ResponseCmd::action(int argc, char **argv)
else
{
// polarizability calculation
// response amplitude nitscf [nite]
assert(argc > 2);
if ( ui->onpe0() )
cout << " ResponseCmd: polarizability with "
......@@ -358,176 +358,9 @@ void ResponseCmd::responseVext(bool rpa, bool ipa, int nitscf, int nite,
int mycol = ctxt->mycol();
Timer tm_comm_drho, tm_write_drho;
if (io == "base64_parallel")
{
#if ! USE_MPI
cout << "cannot use parallel_io when running in serial" << endl;
#endif
// MPI parallel writing with base64 encoding
// all processors on column 1 write collectively
Base64Transcoder xcdr;
int lastproc = nprow - 1;
while ( lastproc >= 0 && ft2.np2_loc(lastproc) == 0 ) lastproc --;
assert(lastproc >= 0);
for (int ispin = 0; ispin < nspin; ispin++)
{
// following part is a modified version of SlaterDet::write in SlaterDet.C
// use group of 3 algorithm to adjust the size of drhor on each processor
tm_comm_drho.start();
int ndiff;
assert(drho_r[ispin].size()==n012loc);
if ( myrow == 0 )
{
ndiff = n012loc % 3;
ctxt->ibcast_send('c',1,1,&ndiff,1);
}
else
{
ctxt->ibcast_recv('c',1,1,&ndiff,1,0,mycol);
}
// send/receive the head/tail of drhor to neighbor
int nsend_left=0, nsend_right=0, nrecv_left=0, nrecv_right=0;
if ( myrow % 3 == 0 )
{
if ( myrow < lastproc )
nsend_right = ndiff;
}
else if ( myrow % 3 == 1 )
{
if ( myrow <= lastproc )
nrecv_left = ndiff;
if ( myrow <= lastproc-1 )
nrecv_right = ndiff;
}
else if ( myrow % 3 == 2 )
{
if ( myrow <= lastproc && myrow > 0 )
nsend_left = ndiff;
}
double rbuf_left[2], rbuf_right[2], sbuf_left[2], sbuf_right[2];
int tmpr_size = n012loc;
if ( nsend_left > 0 )
{
for ( int i = 0; i < ndiff; i++ )
sbuf_left[i] = drho_r[ispin][i];
ctxt->dsend(ndiff,1,sbuf_left,ndiff,myrow-1,mycol);
tmpr_size -= ndiff;
}
if ( nsend_right > 0 )
{
for ( int i = 0; i < ndiff; i++ )
sbuf_right[i] = drho_r[ispin][n012loc-ndiff+i];
ctxt->dsend(ndiff,1,sbuf_right,ndiff,myrow+1,mycol);
tmpr_size -= ndiff;
}
if ( nrecv_left > 0 )
{
ctxt->drecv(ndiff,1,rbuf_left,ndiff,myrow-1,mycol);
tmpr_size += ndiff;
}
if ( nrecv_right > 0 )
{
ctxt->drecv(ndiff,1,rbuf_right,ndiff,myrow+1,mycol);
tmpr_size += ndiff;
}
// at this point, all communications has been finished
// now construct tmpr to store drhor, tmpr has
// size divisible by 3 (except last processor)
if ( myrow < lastproc ) assert(tmpr_size%3 == 0);
vector<double> tmpr(tmpr_size);
if ( nrecv_left > 0 || nrecv_right > 0 )
{
int index = 0;
if ( nrecv_left > 0 )
{
for ( int i = 0; i < ndiff; i++ )
tmpr[index++] = rbuf_left[i];
}
for ( int i = 0; i < n012loc; i++ )
tmpr[index++] = drho_r[ispin][i];
if ( nrecv_right > 0 )
{
for ( int i = 0; i < ndiff; i++ )
tmpr[index++] = rbuf_right[i];
}
assert(index==tmpr_size);
}
else if ( nsend_left > 0 || nsend_right > 0 )
{
int index = 0;
int istart = (nsend_left > 0) ? ndiff : 0;
int iend = (nsend_right > 0) ? n012loc - ndiff : n012loc;
for ( int i = istart; i < iend; i++ )
tmpr[index++] = drho_r[ispin][i];
assert(index==tmpr_size);
}
else
{
for (int i = 0; i < n012loc; i++)
tmpr[i] = drho_r[ispin][i];
assert(tmpr_size==n012loc);
}
// now transform drhor (as stored in tmpr) as base 64 encoding
int nbytes = tmpr.size() * sizeof(double);
int nchars = xcdr.nchars(nbytes);
char* wbuf = new char[nchars];
xcdr.encode(nbytes, (byte *) &tmpr[0], wbuf);
// column 0 write tmpr into file
string filename;
if (nspin == 1)
filename = s->vext->filename() + ".response";
else
filename = s->vext->filename() + ".response."
+ ((ispin == 0) ? "spin0" : "spin1");
tm_comm_drho.stop();
if ( mycol == 0 )
{
// compute offset
int offset;
MPI_Scan(&nchars, &offset, 1, MPI_INT, MPI_SUM, basis.comm());
offset -= nchars;
// write file
tm_write_drho.start();
MPI_File fh;
MPI_Info info;
MPI_Info_create(&info);
int err;
err = MPI_File_open(basis.comm(), (char*) filename.c_str(),
MPI_MODE_WRONLY | MPI_MODE_CREATE, info, &fh);
if (err != 0)
cout << myrow << ": error in MPI_File_open: " << err << endl;
MPI_File_set_size(fh, 0);
MPI_Status status;
err = MPI_File_write_at_all(fh,offset,(void*) &wbuf[0],nchars,
MPI_CHAR,&status);
if ( err != 0 )
cout << myrow << ": error in MPI_File_write_at_all: err="
<< err << endl;
err = MPI_File_close(&fh);
if ( err != 0 )
cout << myrow << ": error in MPI_File_close: " << err << endl;
tm_write_drho.stop();
}
delete [] wbuf;
} // for ispin
}
else if (io == "base64_serial" or io == "cube")
if (io == "xml" or io == "cube")
{
// serial write with base64 encoding or cube format
// serial write xml file or cube file
// processors on row 1 collect drho from other rows
// processor at row 1, column 1 write drho to disk
vector<double> drho_r_gathered;
......@@ -560,17 +393,35 @@ void ResponseCmd::responseVext(bool rpa, bool ipa, int nitscf, int nite,
else
filename = s->vext->filename() + ".response."
+ ((ispin == 0) ? "spin0" : "spin1");
if (io == "base64_serial")
if (io == "xml")
{
Base64Transcoder xcdr;
#if PLT_BIG_ENDIAN
xcdr.byteswap_double(drho_r_gathered.size(),&drho_r_gathered[0]);
#endif
// transform drhor (stored in drho_r_gathered) to base64 encoding
int nbytes = drho_r_gathered.size() * sizeof(double);
int nchars = xcdr.nchars(nbytes);
char *wbuf = new char[nchars];
xcdr.encode(nbytes, (byte *) &drho_r_gathered[0], wbuf);
ofstream os(filename.c_str(), ios::out | ios::binary);
os.write(wbuf, nchars);
ofstream os(filename.c_str(), ios::out );
os << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
os << "<function3d name=\"delta_rho\">" << endl;
D3vector a0 = s->atoms.cell().a(0);
D3vector a1 = s->atoms.cell().a(1);
D3vector a2 = s->atoms.cell().a(2);
os << "<domain a=\"" << a0 << "\"" << endl;
os << " b=\"" << a1 << "\"" << endl;
os << " c=\"" << a2 << "\"/>" << endl;
os << "<grid nx=\"" << n0 << "\" ny=\"" << n1
<< "\" nz=\"" << n2 << "\"/>" << endl;
os << "<grid_function type=\"double\" nx=\"" << n0
<< "\" ny=\"" << n1 << "\" nz=\"" << n2
<< "\" encoding=\"base64\">" << endl;
xcdr.print(nchars,wbuf,os);
os << "</grid_function>" << endl;
os << "</function3d>" << endl;
os.close();
delete [] wbuf;
}
......
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