Homec4science
Diffusion lammps (master) Browse

lammps/src/USER-PHONON/fix_phonon.cpp
c612b724dfc3master

Last Change
Raw File

fix_phonon.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Ling-Ti Kong
Contact:
School of Materials Science and Engineering,
Shanghai Jiao Tong University,
800 Dongchuan Road, Minhang,
Shanghai 200240, CHINA
konglt@sjtu.edu.cn; konglt@gmail.com
------------------------------------------------------------------------- */
#include "string.h"
#include "fix_phonon.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "fft3d_wrap.h"
#include "force.h"
#include "group.h"
#include "lattice.h"
#include "modify.h"
#include "update.h"
#include "citeme.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define MAXLINE 256
static const char cite_fix_phonon[] =
"fix phonon command:\n\n"
"@Article{Kong11,\n"
" author = {L. T. Kong},\n"
" title = {Phonon dispersion measured directly from molecular dynamics simulations},\n"
" journal = {Comp.~Phys.~Comm.},\n"
" year = 2011,\n"
" volume = 182,\n"
" pages = {2201--2207}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_phonon);
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
if (narg < 8) error->all(FLERR,"Illegal fix phonon command: number of arguments < 8");
nevery = force->inumeric(FLERR,arg[3]); // Calculate this fix every n steps!
if (nevery < 1) error->all(FLERR,"Illegal fix phonon command");
nfreq = force->inumeric(FLERR,arg[4]); // frequency to output result
if (nfreq < 1) error->all(FLERR,"Illegal fix phonon command");
waitsteps = ATOBIGINT(arg[5]); // Wait this many timesteps before actually measuring
if (waitsteps < 0) error->all(FLERR,"Illegal fix phonon command: waitsteps < 0 !");
int n = strlen(arg[6]) + 1; // map file
mapfile = new char[n];
strcpy(mapfile, arg[6]);
n = strlen(arg[7]) + 1; // prefix of output
prefix = new char[n];
strcpy(prefix, arg[7]);
logfile = new char[n+4];
sprintf(logfile,"%s.log",prefix);
int sdim = 4;
int iarg = 8;
nasr = 20;
// other command line options
while (iarg < narg){
if (strcmp(arg[iarg],"sysdim") == 0){
if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options.");
sdim = force->inumeric(FLERR,arg[iarg]);
} else if (strcmp(arg[iarg],"nasr") == 0){
if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options.");
nasr = force->inumeric(FLERR,arg[iarg]);
} else {
error->all(FLERR,"Illegal fix phonon command: unknown option read!");
}
iarg++;
}
// get the dimension of the simulation; 1D is possible by specifying the option of "sysdim 1"
sysdim = domain->dimension;
if (sdim < sysdim) sysdim = sdim;
nasr = MAX(0,nasr);
// get the total number of atoms in group
nGFatoms = static_cast<int>(group->count(igroup));
if (nGFatoms<1) error->all(FLERR,"No atom found for fix phonon!");
// MPI gatherv related variables
recvcnts = new int[nprocs];
displs = new int[nprocs];
// mapping index
tag2surf.clear(); // clear map info
surf2tag.clear();
// get the mapping between lattice indices and atom IDs
readmap(); delete []mapfile;
if (nucell == 1) nasr = MIN(1,nasr);
// get the mass matrix for dynamic matrix
getmass();
// create FFT and allocate memory for FFT
nxlo = 0;
int *nx_loc = new int [nprocs];
for (int i=0; i<nprocs;i++){
nx_loc[i] = nx/nprocs;
if (i < nx%nprocs) nx_loc[i]++;
}
for (int i=0; i<me; i++) nxlo += nx_loc[i];
nxhi = nxlo + nx_loc[me] - 1;
mynpt = nx_loc[me]*ny*nz;
mynq = mynpt;
fft_dim = nucell*sysdim;
fft_dim2 = fft_dim*fft_dim;
fft_nsend = mynpt*fft_dim;
fft_cnts = new int[nprocs];
fft_disp = new int[nprocs];
fft_disp[0] = 0;
for (int i=0; i<nprocs; i++) fft_cnts[i] = nx_loc[i]*ny*nz*fft_dim;
for (int i=1; i<nprocs; i++) fft_disp[i] = fft_disp[i-1]+fft_cnts[i-1];
delete []nx_loc;
fft = new FFT3d(lmp,world,nz,ny,nx,0,nz-1,0,ny-1,nxlo,nxhi,0,nz-1,0,ny-1,nxlo,nxhi,0,0,&mysize);
fft_data = (double *) memory->smalloc(MAX(1,mynq)*2*sizeof(double),"fix_phonon:fft_data");
// allocate variables; MAX(1,... is used because NULL buffer will result in error for MPI
RIloc = memory->create(RIloc,nGFatoms,(sysdim+1),"fix_phonon:RIloc");
RIall = memory->create(RIall,nGFatoms,(sysdim+1),"fix_phonon:RIall");
Rsort = memory->create(Rsort,nGFatoms, sysdim, "fix_phonon:Rsort");
Rnow = memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow");
Rsum = memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum");
basis = memory->create(basis,nucell, sysdim, "fix_phonon:basis");
// because of hermit, only nearly half of q points are stored
Rqnow = memory->create(Rqnow,MAX(1,mynq),fft_dim, "fix_phonon:Rqnow");
Rqsum = memory->create(Rqsum,MAX(1,mynq),fft_dim2,"fix_phonon:Rqsum");
Phi_q = memory->create(Phi_q,MAX(1,mynq),fft_dim2,"fix_phonon:Phi_q");
if (me == 0) // variable to collect all local Phi to root
Phi_all = memory->create(Phi_all,ntotal,fft_dim2,"fix_phonon:Phi_all");
else
Phi_all = memory->create(Phi_all,1,1,"fix_phonon:Phi_all");
// output some information on the system to log file
if (me == 0){
flog = fopen(logfile, "w");
if (flog == NULL) {
char str[MAXLINE];
sprintf(str,"Can not open output file %s",logfile);
error->one(FLERR,str);
}
for (int i=0; i<60; i++) fprintf(flog,"#"); fprintf(flog,"\n");
fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]);
fprintf(flog,"# total number of atoms in the group : %d\n", nGFatoms);
fprintf(flog,"# dimension of the system : %d D\n", sysdim);
fprintf(flog,"# number of atoms per unit cell : %d\n", nucell);
fprintf(flog,"# dimension of the FFT mesh : %d x %d x %d\n", nx, ny, nz);
fprintf(flog,"# number of wait steps before measurement : " BIGINT_FORMAT "\n", waitsteps);
fprintf(flog,"# frequency of GFC measurement : %d\n", nevery);
fprintf(flog,"# output result after this many measurement: %d\n", nfreq);
fprintf(flog,"# number of processors used by this run : %d\n", nprocs);
for (int i=0; i<60; i++) fprintf(flog,"#"); fprintf(flog,"\n");
fprintf(flog,"# mapping information between lattice index and atom id\n");
fprintf(flog,"# nx ny nz nucell\n");
fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell);
fprintf(flog,"# l1 l2 l3 k atom_id\n");
int ix, iy, iz, iu;
for (idx =0; idx<nGFatoms; idx++){
itag = surf2tag[idx];
iu = idx%nucell;
iz = (idx/nucell)%nz;
iy = (idx/(nucell*nz))%ny;
ix = (idx/(nucell*nz*ny))%nx;
fprintf(flog,"%d %d %d %d %d\n", ix, iy, iz, iu, itag);
}
for (int i=0; i<60; i++) fprintf(flog,"#"); fprintf(flog,"\n");
fflush(flog);
}
surf2tag.clear();
// default temperature is from thermo
TempSum = new double[sysdim];
id_temp = new char[12];
strcpy(id_temp,"thermo_temp");
int icompute = modify->find_compute(id_temp);
temperature = modify->compute[icompute];
inv_nTemp = 1.0/group->count(temperature->igroup);
} // end of constructor
/* ---------------------------------------------------------------------- */
void FixPhonon::post_run()
{
// compute and output final GFC results
if (ifreq > 0 && ifreq != nfreq) postprocess();
if (me == 0) fclose(flog);
}
/* ---------------------------------------------------------------------- */
FixPhonon::~FixPhonon()
{
// delete locally stored array
memory->destroy(RIloc);
memory->destroy(RIall);
memory->destroy(Rsort);
memory->destroy(Rnow);
memory->destroy(Rsum);
memory->destroy(basis);
memory->destroy(Rqnow);
memory->destroy(Rqsum);
memory->destroy(Phi_q);
memory->destroy(Phi_all);
delete []recvcnts;
delete []displs;
delete []prefix;
delete []logfile;
delete []fft_cnts;
delete []fft_disp;
delete []id_temp;
delete []TempSum;
delete []M_inv_sqrt;
delete []basetype;
// destroy FFT
delete fft;
memory->sfree(fft_data);
// clear map info
tag2surf.clear();
surf2tag.clear();
}
/* ---------------------------------------------------------------------- */
int FixPhonon::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPhonon::init()
{
// warn if more than one fix-phonon
int count = 0;
for (int i=0;i<modify->nfix;i++) if (strcmp(modify->fix[i]->style,"gfc") == 0) count++;
if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed.
}
/* ---------------------------------------------------------------------- */
void FixPhonon::setup(int flag)
{
// initialize accumulating variables
for (int i=0; i<sysdim; i++) TempSum[i] = 0.;
for (int i=0; i<mynpt; i++){
for (int j=0; j<fft_dim; j++) Rsum[i][j] = 0.;
}
for (int i=0; i<mynq; i++){
for (int j=0; j<fft_dim2; j++) Rqsum[i][j] = std::complex<double> (0.,0.);
}
for (int i=0; i<6; i++) hsum[i] = 0.;
for (int i=0; i<nucell; i++){
for (int j=0; j<sysdim; j++) basis[i][j] = 0.;
}
prev_nstep = update->ntimestep;
GFcounter = 0;
ifreq = 0;
}
/* ---------------------------------------------------------------------- */
void FixPhonon::end_of_step()
{
if ( (update->ntimestep-prev_nstep) <= waitsteps) return;
double **x = atom->x;
int *mask = atom->mask;
int *tag = atom->tag;
int *image = atom->image;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double *h = domain->h;
double xbox, ybox, zbox;
int i,idim,jdim,ndim;
double xcur[3];
// to get the current temperature
if (!(temperature->invoked_flag & INVOKED_VECTOR)) temperature->compute_vector();
for (idim=0; idim<sysdim; idim++) TempSum[idim] += temperature->vector[idim];
// evaluate R(r) on local proc
nfind = 0;
if (domain->triclinic == 0) { // for orthogonal lattice
for (i=0; i<nlocal; i++){
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
xcur[0] = x[i][0] + xprd*xbox;
xcur[1] = x[i][1] + yprd*ybox;
xcur[2] = x[i][2] + zprd*zbox;
for (idim=0; idim<sysdim; idim++) RIloc[nfind][idim] = xcur[idim];
RIloc[nfind++][sysdim] = idx;
}
}
}else{ // for non-orthogonal lattice
for (i=0; i<nlocal; i++){
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
xcur[0] = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
xcur[1] = x[i][1] + h[1]*ybox + h[3]*zbox;
xcur[2] = x[i][2] + h[2]*zbox;
for (idim=0; idim<sysdim; idim++) RIloc[nfind][idim] = xcur[idim];
RIloc[nfind++][sysdim] = idx;
}
}
}
// gather R(r) on local proc, then sort and redistribute to all procs for FFT
nfind *= (sysdim+1);
displs[0] = 0;
for (i=0; i<nprocs;i++) recvcnts[i] = 0;
MPI_Gather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,0,world);
for (i=1; i<nprocs; i++) displs[i] = displs[i-1] + recvcnts[i-1];
MPI_Gatherv(RIloc[0],nfind,MPI_DOUBLE,RIall[0],recvcnts,displs,MPI_DOUBLE,0,world);
if (me == 0){
for (i=0; i<nGFatoms; i++){
idx = static_cast<int>(RIall[i][sysdim]);
for (idim=0; idim<sysdim; idim++) Rsort[idx][idim] = RIall[i][idim];
}
}
MPI_Scatterv(Rsort[0],fft_cnts,fft_disp, MPI_DOUBLE, Rnow[0], fft_nsend, MPI_DOUBLE,0,world);
// get Rsum
for (idx=0; idx<mynpt; idx++){
for (idim=0; idim<fft_dim; idim++) Rsum[idx][idim] += Rnow[idx][idim];
}
// FFT R(r) to get R(q)
for (idim=0; idim<fft_dim; idim++){
int m=0;
for (idx=0; idx<mynpt; idx++){
fft_data[m++] = Rnow[idx][idim];
fft_data[m++] = 0.;
}
fft->compute(fft_data, fft_data, -1);
m = 0;
for (idq=0; idq<mynq; idq++){
Rqnow[idq][idim] = std::complex<double>(fft_data[m], fft_data[m+1]);
m += 2;
}
}
// to get sum(R(q).R(q)*)
for (idq=0; idq<mynq; idq++){
ndim = 0;
for (idim=0; idim<fft_dim; idim++){
for (jdim=0; jdim<fft_dim; jdim++){
Rqsum[idq][ndim++] += Rqnow[idq][idim]*conj(Rqnow[idq][jdim]);
}
}
}
// get basis info
if (fft_dim > sysdim){
for (idx=0; idx<mynpt; idx++){
ndim = sysdim;
for (i=1; i<nucell; i++){
double dist2orig[3];
for (idim=0; idim<sysdim; idim++) dist2orig[idim] = Rnow[idx][ndim++] - Rnow[idx][idim];
domain->minimum_image(dist2orig);
for (idim=0; idim<sysdim; idim++) basis[i][idim] += dist2orig[idim];
}
}
}
// get lattice vector info
for (int i=0; i<6; i++) hsum[i] += h[i];
// increment counter
GFcounter++;
// compute and output Phi_q after every nfreq evaluations
if (++ifreq == nfreq) postprocess();
} // end of end_of_step()
/* ---------------------------------------------------------------------- */
double FixPhonon::memory_usage()
{
double bytes = sizeof(double)*2*mynq
+ sizeof(std::map<int,int>)*2*nGFatoms
+ sizeof(double)*(nGFatoms*(3*sysdim+2)+mynpt*fft_dim*2)
+ sizeof(std::complex<double>)*MAX(1,mynq)*fft_dim *(1+2*fft_dim)
+ sizeof(std::complex<double>)*ntotal*fft_dim2
+ sizeof(int) * nprocs * 4;
return bytes;
}
/* ---------------------------------------------------------------------- */
int FixPhonon::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify temp ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,"Fix_modify temp ID does not compute temperature");
inv_nTemp = 1.0/group->count(temperature->igroup);
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
* private method, to get the mass matrix for dynamic matrix
* --------------------------------------------------------------------*/
void FixPhonon::getmass()
{
int nlocal = atom->nlocal;
int *mask = atom->mask;
int *tag = atom->tag;
int *type = atom->type;
double *rmass = atom->rmass;
double *mass = atom->mass;
double *mass_one, *mass_all;
double *type_one, *type_all;
mass_one = new double[nucell];
mass_all = new double[nucell];
type_one = new double[nucell];
type_all = new double[nucell];
for (int i=0; i<nucell; i++) {
mass_one[i] = 0.;
type_one[i] = 0.;
}
if (rmass){
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
int iu = idx%nucell;
mass_one[iu] += rmass[i];
type_one[iu] += double(type[i]);
}
}
} else {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
int iu = idx%nucell;
mass_one[iu] += mass[type[i]];
type_one[iu] += double(type[i]);
}
}
}
MPI_Allreduce(mass_one,mass_all,nucell,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(type_one,type_all,nucell,MPI_DOUBLE,MPI_SUM,world);
M_inv_sqrt = new double[nucell];
basetype = new int[nucell];
double inv_total = 1./double(ntotal);
for (int i=0; i<nucell; i++){
mass_all[i] *= inv_total;
M_inv_sqrt[i] = sqrt(1./mass_all[i]);
basetype[i] = int(type_all[i]*inv_total);
}
delete []mass_one;
delete []mass_all;
delete []type_one;
delete []type_all;
return;
}
/* ----------------------------------------------------------------------
* private method, to read the mapping info from file
* --------------------------------------------------------------------*/
void FixPhonon::readmap()
{
int info = 0;
char strtmp[MAXLINE];
FILE *fp;
fp = fopen(mapfile, "r");
if (fp == NULL){
char str[MAXLINE];
sprintf(str,"Cannot open input map file %s", mapfile);
error->all(FLERR,str);
}
if (fgets(strtmp,MAXLINE,fp) == NULL)
error->all(FLERR,"Error while reading header of mapping file!");
sscanf(strtmp,"%d %d %d %d", &nx, &ny, &nz, &nucell);
ntotal = nx*ny*nz;
if (ntotal*nucell != nGFatoms) error->all(FLERR,"FFT mesh and number of atoms in group mismatch!");
if (fgets(strtmp,MAXLINE,fp) == NULL) // second line of mapfile is comment
error->all(FLERR,"Error while reading comment of mapping file!");
int ix, iy, iz, iu;
for (int i=0; i<nGFatoms; i++){
if (fgets(strtmp,MAXLINE,fp) == NULL) {info = 1; break;}
sscanf(strtmp,"%d %d %d %d %d", &ix, &iy, &iz, &iu, &itag); // the remaining lines carry the mapping info
if (ix<0 || ix>=nx || iy<0 || iy>=ny || iz<0 || iz>=nz|| iu<0 || iu>=nucell) {info = 2; break;} // check if index is in correct range
if (itag<1 || itag>static_cast<int>(atom->natoms)) {info = 3; break;} // 1 <= itag <= natoms
idx = ((ix*ny+iy)*nz+iz)*nucell+iu;
tag2surf[itag] = idx;
surf2tag[idx] = itag;
}
fclose(fp);
if (tag2surf.size() != surf2tag.size() || tag2surf.size() != static_cast<std::size_t>(nGFatoms) )
error->all(FLERR,"The mapping is incomplete!");
if (info) error->all(FLERR,"Error while reading mapping file!");
// check the correctness of mapping
int *mask = atom->mask;
int *tag = atom->tag;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
if (itag != surf2tag[idx]) error->one(FLERR,"The mapping info read is incorrect!");
}
}
return;
}
/* ----------------------------------------------------------------------
* private method, to output the force constant matrix
* --------------------------------------------------------------------*/
void FixPhonon::postprocess( )
{
if (GFcounter<1) return;
ifreq =0;
int idim, jdim, ndim;
double invGFcounter = 1.0 /double(GFcounter);
// to get <Rq.Rq*>
for (idq=0; idq<mynq; idq++){
for (idim=0; idim<fft_dim2; idim++) Phi_q[idq][idim] = Rqsum[idq][idim]*invGFcounter;
}
// to get <R>
for (idx=0; idx<mynpt; idx++){
for (idim=0; idim<fft_dim; idim++) Rnow[idx][idim] = Rsum[idx][idim] * invGFcounter;
}
// to get <R>q
for (idim=0; idim<fft_dim; idim++){
int m = 0;
for (idx=0; idx<mynpt; idx++){
fft_data[m++] = Rnow[idx][idim];
fft_data[m++] = 0.;
}
fft->compute(fft_data,fft_data,-1);
m = 0;
for (idq=0; idq<mynq; idq++){
Rqnow[idq][idim] = std::complex<double>(fft_data[m], fft_data[m+1]);
m += 2;
}
}
// to get G(q) = <Rq.Rq*> - <R>q.<R*>q
for (idq=0; idq<mynq; idq++){
ndim = 0;
for (idim=0; idim<fft_dim; idim++){
for (jdim=0; jdim<fft_dim; jdim++) Phi_q[idq][ndim++] -= Rqnow[idq][idim]*conj(Rqnow[idq][jdim]);
}
}
// to get Phi = KT.G^-1; normalization of FFTW data is done here
double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.;
double TempFac = invGFcounter*inv_nTemp;
double NormFac = TempFac*double(ntotal);
for (idim=0; idim<sysdim; idim++){
kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac);
TempAve += TempSum[idim]*TempFac;
}
TempAve /= sysdim*boltz;
for (idq=0; idq<mynq; idq++){
GaussJordan(fft_dim, Phi_q[idq]);
ndim =0;
for (idim=0; idim<fft_dim; idim++){
for (jdim=0; jdim<fft_dim; jdim++){
Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim];
}
}
}
// to collect all local Phi_q to root
displs[0]=0;
for (int i=0; i<nprocs; i++) recvcnts[i] = fft_cnts[i]*fft_dim*2;
for (int i=1; i<nprocs; i++) displs[i] = displs[i-1] + recvcnts[i-1];
MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world);
// to collect all basis info and averaged it on root
double basis_root[fft_dim];
if (fft_dim > sysdim)
MPI_Reduce (&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world);
if (me == 0){ // output dynamic matrix by root
// get basis info
for (idim=0; idim<sysdim; idim++) basis_root[idim] = 0.;
for (idim=sysdim; idim<fft_dim; idim++) basis_root[idim] /= double(ntotal)*double(GFcounter);
// get unit cell base vector info; might be incorrect if MD pbc and FixPhonon pbc mismatch.
double basevec[9];
basevec[1] = basevec[2] = basevec[5] = 0.;
basevec[0] = hsum[0] * invGFcounter / double(nx);
basevec[4] = hsum[1] * invGFcounter / double(ny);
basevec[8] = hsum[2] * invGFcounter / double(nz);
basevec[7] = hsum[3] * invGFcounter / double(nz);
basevec[6] = hsum[4] * invGFcounter / double(nz);
basevec[3] = hsum[5] * invGFcounter / double(ny);
// write binary file, in fact, it is the force constants matrix that is written
// Enforcement of ASR and the conversion of dynamical matrix is done in the postprocessing code
char fname[MAXLINE];
FILE *GFC_bin;
sprintf(fname,"%s.bin." BIGINT_FORMAT,prefix,update->ntimestep);
GFC_bin = fopen(fname,"wb");
fwrite(&sysdim, sizeof(int), 1, GFC_bin);
fwrite(&nx, sizeof(int), 1, GFC_bin);
fwrite(&ny, sizeof(int), 1, GFC_bin);
fwrite(&nz, sizeof(int), 1, GFC_bin);
fwrite(&nucell, sizeof(int), 1, GFC_bin);
fwrite(&boltz, sizeof(double), 1, GFC_bin);
fwrite(Phi_all[0],sizeof(double),ntotal*fft_dim2*2,GFC_bin);
fwrite(&TempAve, sizeof(double),1, GFC_bin);
fwrite(&basevec[0], sizeof(double),9, GFC_bin);
fwrite(&basis_root[0],sizeof(double),fft_dim,GFC_bin);
fwrite(basetype, sizeof(int), nucell, GFC_bin);
fwrite(M_inv_sqrt, sizeof(double),nucell, GFC_bin);
fclose(GFC_bin);
// write log file, here however, it is the dynamical matrix that is written
for (int i=0; i<60; i++) fprintf(flog,"#"); fprintf(flog,"\n");
fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep);
fprintf(flog, "# Total number of measurements : %d\n", GFcounter);
fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve);
fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz);
fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]);
fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]);
fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]);
for (int i=0; i<60; i++) fprintf(flog,"#"); fprintf(flog,"\n");
fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n");
EnforceASR();
// to get D = 1/M x Phi
for (idq=0; idq<ntotal; idq++){
ndim =0;
for (idim=0; idim<fft_dim; idim++){
for (jdim=0; jdim<fft_dim; jdim++){
Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
}
}
}
idq =0;
for (int ix=0; ix<nx; ix++){
double qx = double(ix)/double(nx);
for (int iy=0; iy<ny; iy++){
double qy = double(iy)/double(ny);
for (int iz=0; iz<nz; iz++){
double qz = double(iz)/double(nz);
fprintf(flog,"%lg %lg %lg", qx, qy, qz);
for (idim=0; idim<fft_dim2; idim++) fprintf(flog, " %lg %lg", real(Phi_all[idq][idim]), imag(Phi_all[idq][idim]));
fprintf(flog, "\n");
idq++;
}
}
}
fflush(flog);
}
} // end of postprocess
/* ----------------------------------------------------------------------
* private method, to get the inverse of a complex matrix by means of
* Gaussian-Jordan Elimination with full pivoting; square matrix required.
*
* Adapted from the Numerical Recipes in Fortran.
* --------------------------------------------------------------------*/
void FixPhonon::GaussJordan(int n, std::complex<double> *Mat)
{
int i,icol,irow,j,k,l,ll,idr,idc;
int *indxc,*indxr,*ipiv;
double big, nmjk;
std::complex<double> dum, pivinv;
indxc = new int[n];
indxr = new int[n];
ipiv = new int[n];
for (i=0; i<n; i++) ipiv[i] = 0;
for (i=0; i<n; i++){
big = 0.;
for (j=0; j<n; j++){
if (ipiv[j] != 1){
for (k=0; k<n; k++){
if (ipiv[k] == 0){
idr = j*n+k;
nmjk = norm(Mat[idr]);
if (nmjk >= big){
big = nmjk;
irow = j;
icol = k;
}
}else if (ipiv[k]>1){
error->one(FLERR,"Singular matrix in complex GaussJordan!");
}
}
}
}
ipiv[icol] += 1;
if (irow != icol){
for (l=0; l<n; l++){
idr = irow*n+l;
idc = icol*n+l;
dum = Mat[idr];
Mat[idr] = Mat[idc];
Mat[idc] = dum;
}
}
indxr[i] = irow;
indxc[i] = icol;
idr = icol*n+icol;
if (Mat[idr] == std::complex<double>(0.,0.)) error->one(FLERR,"Singular matrix in complex GaussJordan!");
pivinv = 1./ Mat[idr];
Mat[idr] = std::complex<double>(1.,0.);
idr = icol*n;
for (l=0; l<n; l++) Mat[idr+l] *= pivinv;
for (ll=0; ll<n; ll++){
if (ll != icol){
idc = ll*n+icol;
dum = Mat[idc];
Mat[idc] = 0.;
idc -= icol;
for (l=0; l<n; l++) Mat[idc+l] -= Mat[idr+l]*dum;
}
}
}
for (l=n-1; l>=0; l--){
int rl = indxr[l];
int cl = indxc[l];
if (rl != cl){
for (k=0; k<n; k++){
idr = k*n+rl;
idc = k*n+cl;
dum = Mat[idr];
Mat[idr] = Mat[idc];
Mat[idc] = dum;
}
}
}
delete []indxr;
delete []indxc;
delete []ipiv;
return;
}
/* ----------------------------------------------------------------------
* private method, to apply the acoustic sum rule on force constant matrix
* at gamma point. Should be executed on root only.
* --------------------------------------------------------------------*/
void FixPhonon::EnforceASR()
{
if (nasr < 1) return;
for (int iit=0; iit<nasr; iit++){
// simple ASR; the resultant matrix might not be symmetric
for (int a=0; a<sysdim; a++)
for (int b=0; b<sysdim; b++){
for (int k=0; k<nucell; k++){
double sum = 0.;
for (int kp=0; kp<nucell; kp++){
int idx = (k*sysdim+a)*fft_dim+kp*sysdim+b;
sum += real(Phi_all[0][idx]);
}
sum /= double(nucell);
for (int kp=0; kp<nucell; kp++){
int idx = (k*sysdim+a)*fft_dim+kp*sysdim+b;
real(Phi_all[0][idx]) -= sum;
}
}
}
// symmetrize
for (int k=0; k<nucell; k++)
for (int kp=k; kp<nucell; kp++){
double csum = 0.;
for (int a=0; a<sysdim; a++)
for (int b=0; b<sysdim; b++){
int idx = (k*sysdim+a)*fft_dim+kp*sysdim+b;
int jdx = (kp*sysdim+b)*fft_dim+k*sysdim+a;
csum = (real(Phi_all[0][idx])+real(Phi_all[0][jdx]))*0.5;
real(Phi_all[0][idx]) = real(Phi_all[0][jdx]) = csum;
}
}
}
// symmetric ASR
for (int a=0; a<sysdim; a++)
for (int b=0; b<sysdim; b++){
for (int k=0; k<nucell; k++){
double sum = 0.;
for (int kp=0; kp<nucell; kp++){
int idx = (k*sysdim+a)*fft_dim+kp*sysdim+b;
sum += real(Phi_all[0][idx]);
}
sum /= double(nucell-k);
for (int kp=k; kp<nucell; kp++){
int idx = (k*sysdim+a)*fft_dim+kp*sysdim+b;
int jdx = (kp*sysdim+b)*fft_dim+k*sysdim+a;
real(Phi_all[0][idx]) -= sum;
real(Phi_all[0][jdx]) = real(Phi_all[0][idx]);
}
}
}
}
/* --------------------------------------------------------------------*/
c4science · Help