git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14372 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-12-15 22:30:45 +00:00
parent b5a1ba9bfa
commit 0939eb1ee9
23 changed files with 88 additions and 90 deletions

View File

@ -25,7 +25,7 @@
struct lmp_float3 {
float x,y,z;
KOKKOS_INLINE_FUNCTION
lmp_float3():x(0.0f),z(0.0f),y(0.0f) {}
lmp_float3():x(0.0f),y(0.0f),z(0.0f) {}
KOKKOS_INLINE_FUNCTION
void operator += (const lmp_float3& tmp) {
@ -56,7 +56,7 @@
struct lmp_double3 {
double x,y,z;
KOKKOS_INLINE_FUNCTION
lmp_double3():x(0.0),z(0.0),y(0.0) {}
lmp_double3():x(0.0),y(0.0),z(0.0) {}
KOKKOS_INLINE_FUNCTION
void operator += (const lmp_double3& tmp) {
@ -744,5 +744,9 @@ void memset_kokkos (ViewType &view) {
ViewType::execution_space::fence();
}
#if defined(KOKKOS_HAVE_CXX11)
#undef ISFINITE
#define ISFINITE(x) std::isfinite(x)
#endif
#endif

View File

@ -103,7 +103,6 @@ void PairCoulDSFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
q = atomKK->k_q.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
@ -250,7 +249,6 @@ void PairCoulDSFKokkos<DeviceType>::operator()(TagPairCoulDSFKernelA<NEIGHFLAG,N
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const F_FLOAT qtmp = q[i];
const int itype = type(i);
if (eflag) {
const F_FLOAT e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
@ -275,7 +273,6 @@ void PairCoulDSFKokkos<DeviceType>::operator()(TagPairCoulDSFKernelA<NEIGHFLAG,N
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
const int jtype = type(j);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cut_coulsq) {
@ -432,4 +429,4 @@ int PairCoulDSFKokkos<DeviceType>::sbmask(const int& j) const {
template class PairCoulDSFKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairCoulDSFKokkos<LMPHostType>;
#endif
#endif

View File

@ -66,7 +66,6 @@ class PairCoulDSFKokkos : public PairCoulDSF {
typename ArrayTypes<DeviceType>::t_x_array_randomread x;
typename ArrayTypes<DeviceType>::t_f_array f;
typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;

View File

@ -104,7 +104,6 @@ void PairCoulWolfKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
q = atomKK->k_q.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
@ -251,7 +250,6 @@ void PairCoulWolfKokkos<DeviceType>::operator()(TagPairCoulWolfKernelA<NEIGHFLAG
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const F_FLOAT qtmp = q[i];
const int itype = type(i);
if (eflag) {
const F_FLOAT qisq = qtmp*qtmp;
@ -277,7 +275,6 @@ void PairCoulWolfKokkos<DeviceType>::operator()(TagPairCoulWolfKernelA<NEIGHFLAG
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
const int jtype = type(j);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cut_coulsq) {
@ -434,4 +431,4 @@ int PairCoulWolfKokkos<DeviceType>::sbmask(const int& j) const {
template class PairCoulWolfKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairCoulWolfKokkos<LMPHostType>;
#endif
#endif

View File

@ -66,7 +66,6 @@ class PairCoulWolfKokkos : public PairCoulWolf {
typename ArrayTypes<DeviceType>::t_x_array_randomread x;
typename ArrayTypes<DeviceType>::t_f_array f;
typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;

View File

@ -1432,7 +1432,7 @@ void MSM::particle_map()
int flag = 0;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int i = 0; i < nlocal; i++) {

View File

@ -314,7 +314,7 @@ void MSMCG::particle_map()
int flag = 0;
int i;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int j = 0; j < num_charged; j++) {

View File

@ -1890,7 +1890,7 @@ void PPPM::particle_map()
int flag = 0;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
for (int i = 0; i < nlocal; i++) {

View File

@ -282,7 +282,7 @@ void PPPMCG::particle_map()
double **x = atom->x;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;

View File

@ -4229,7 +4229,7 @@ void PPPMDisp::particle_map(double delx, double dely, double delz,
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;

View File

@ -78,7 +78,7 @@ void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz,
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;

View File

@ -679,7 +679,7 @@ void PPPMStagger::particle_map()
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;

View File

@ -73,7 +73,7 @@ void PPPMTIP4P::particle_map()
double **x = atom->x;
int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int flag = 0;

View File

@ -331,7 +331,7 @@ void MSMCGOMP::particle_map()
int flag = 0;
int i;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
// XXX: O(N). is it worth to add OpenMP here?

View File

@ -361,7 +361,7 @@ void PPPMDispOMP::particle_map(double dxinv, double dyinv,
const int nyhi_out = nyhi_o;
const int nzhi_out = nzhi_o;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions. Simulation unstable.");
int i, flag = 0;

View File

@ -359,7 +359,7 @@ void PPPMDispTIP4POMP::particle_map_c(double dxinv, double dyinv,
const int nyhi_out = nyhi_o;
const int nzhi_out = nzhi_o;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int i, flag = 0;

View File

@ -350,7 +350,7 @@ void PPPMTIP4POMP::particle_map()
const double boxloz = boxlo[2];
const int nlocal = atom->nlocal;
if (!isfinite(boxlo[0]) || !isfinite(boxlo[1]) || !isfinite(boxlo[2]))
if (!ISFINITE(boxlo[0]) || !ISFINITE(boxlo[1]) || !ISFINITE(boxlo[2]))
error->one(FLERR,"Non-numeric box dimensions - simulation unstable");
int i, flag = 0;

View File

@ -15,6 +15,12 @@ There is an auxiliary post-processing tool in tools/phonon that will
compute phonon frequencies and dispersion relations from the dynamical
matrices output by this command.
There is also an alternative code, dump2phonon, available which enables
one to use the functions of fix-phonon by reading in atom-style dump
files of lammps (which can be converted from the trajectories of any
other MD code):
https://github.com/lingtikong/dump2phonon
The person who created this package is Ling-Ti Kong (konglt at
sjtu.edu.cn) at Shanghai Jiao Tong University. Contact him directly
if you have questions.

View File

@ -28,7 +28,6 @@
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "fft3d_wrap.h"
#include "force.h"
#include "group.h"
#include "lattice.h"
@ -64,7 +63,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
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!
@ -85,7 +84,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
strcpy(prefix, arg[7]);
logfile = new char[n+4];
sprintf(logfile,"%s.log",prefix);
int sdim = sysdim = domain->dimension;
int iarg = 8;
nasr = 20;
@ -139,23 +138,23 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nxlo = 0;
int *nx_loc = new int [nprocs];
for (int i = 0; i < nprocs; ++i){
nx_loc[i] = nx/nprocs;
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;
mynpt = nx_loc[me] * ny * nz;
mynq = mynpt;
fft_dim = nucell*sysdim;
fft_dim2 = fft_dim*fft_dim;
fft_nsend = mynpt*fft_dim;
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];
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,0);
@ -165,10 +164,10 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
memory->create(RIloc,ngroup,(sysdim+1),"fix_phonon:RIloc");
memory->create(RIall,ngroup,(sysdim+1),"fix_phonon:RIall");
memory->create(Rsort,ngroup, sysdim, "fix_phonon:Rsort");
memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow");
memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum");
memory->create(basis,nucell, sysdim, "fix_phonon:basis");
// because of hermit, only nearly half of q points are stored
@ -199,7 +198,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
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,"# mapping information between lattice indices 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");
@ -216,14 +215,14 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
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);
inv_nTemp = 1./group->count(temperature->igroup);
} // end of constructor
@ -268,7 +267,7 @@ FixPhonon::~FixPhonon()
// destroy FFT
delete fft;
memory->sfree(fft_data);
// clear map info
tag2surf.clear();
surf2tag.clear();
@ -291,7 +290,7 @@ 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;
for (int i = 0; i < modify->nfix; ++i) if (strcmp(modify->fix[i]->style,"phonon") == 0) ++count;
if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed.
}
@ -313,8 +312,8 @@ void FixPhonon::setup(int flag)
for (int i = 0; i < nucell; ++i)
for (int j = 0; j < sysdim; ++j) basis[i][j] = 0.;
neval = ifreq = 0;
prev_nstep = update->ntimestep;
neval = ifreq = 0;
}
/* ---------------------------------------------------------------------- */
@ -324,12 +323,12 @@ void FixPhonon::end_of_step()
if ( (update->ntimestep-prev_nstep) <= waitsteps) return;
double **x = atom->x;
int *mask = atom->mask;
tagint *tag = atom->tag;
int *mask = atom->mask;
tagint *tag = atom->tag;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double *h = domain->h;
double *h = domain->h;
int i,idim,jdim,ndim;
double xcur[3];
@ -346,9 +345,9 @@ void FixPhonon::end_of_step()
idx = tag2surf[itag];
domain->unmap(x[i], image[i], xcur);
for (idim = 0; idim < sysdim; ++idim) RIloc[nfind][idim] = xcur[idim];
RIloc[nfind++][sysdim] = double(idx);
RIloc[nfind++][sysdim] = static_cast<double>(idx);
}
}
@ -374,17 +373,17 @@ void FixPhonon::end_of_step()
// FFT R(r) to get R(q)
for (idim = 0; idim < fft_dim; ++idim){
int m=0;
int m = 0;
for (idx = 0; idx < mynpt; ++idx){
fft_data[m++] = Rnow[idx][idim];
fft_data[m++] = 0.;
fft_data[m++] = static_cast<FFT_SCALAR>(Rnow[idx][idim]);
fft_data[m++] = static_cast<FFT_SCALAR>(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]);
Rqnow[idq][idim] = std::complex<double>(static_cast<double>(fft_data[m]), static_cast<double>(fft_data[m+1]));
m += 2;
}
}
@ -393,7 +392,7 @@ void FixPhonon::end_of_step()
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]);
for (jdim = 0; jdim < fft_dim; ++jdim) Rqsum[idq][ndim++] += Rqnow[idq][idim] * std::conj(Rqnow[idq][jdim]);
}
// get basis info
@ -540,14 +539,14 @@ void FixPhonon::readmap()
for (int i = 0; i < atom->nlocal; ++i){
if (atom->mask[i] & groupbit) tag_loc[nfind++] = atom->tag[i];
}
// gather IDs on local proc
displs[0] = 0;
for (int i = 0; i < nprocs; ++i) recvcnts[i] = 0;
MPI_Allgather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,world);
for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1];
MPI_Allgatherv(tag_loc,nfind,MPI_LMP_TAGINT,tag_all,recvcnts,displs,MPI_INT,world);
MPI_Allgatherv(tag_loc,nfind,MPI_LMP_TAGINT,tag_all,recvcnts,displs,MPI_LMP_TAGINT,world);
for (int i = 0; i < ngroup; ++i){
itag = tag_all[i];
tag2surf[itag] = i;
@ -567,24 +566,24 @@ void FixPhonon::readmap()
error->all(FLERR,line);
}
if (fgets(line,MAXLINE,fp) == NULL)
if (fgets(line,MAXLINE,fp) == NULL)
error->all(FLERR,"Error while reading header of mapping file!");
nx = force->inumeric(FLERR, strtok(line, " \n\t\r\f"));
ny = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
nz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
nucell = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
ntotal = nx*ny*nz;
if (ntotal*nucell != ngroup)
if (ntotal*nucell != ngroup)
error->all(FLERR,"FFT mesh and number of atoms in group mismatch!");
// second line of mapfile is comment
if (fgets(line,MAXLINE,fp) == NULL)
if (fgets(line,MAXLINE,fp) == NULL)
error->all(FLERR,"Error while reading comment of mapping file!");
int ix, iy, iz, iu;
// the remaining lines carry the mapping info
for (int i = 0; i < ngroup; ++i){
if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;}
if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;}
ix = force->inumeric(FLERR, strtok(line, " \n\t\r\f"));
iy = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
iz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
@ -592,21 +591,21 @@ void FixPhonon::readmap()
itag = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f"));
// check if index is in correct range
if (ix < 0 || ix >= nx || iy < 0 || iy >= ny ||
if (ix < 0 || ix >= nx || iy < 0 || iy >= ny ||
iz < 0 || iz >= nz || iu < 0 || iu >= nucell) {info = 2; break;}
// 1 <= itag <= natoms
if (itag < 1 || itag > static_cast<int>(atom->natoms)) {info = 3; break;}
if (itag < 1 || itag > static_cast<tagint>(atom->natoms)) {info = 3; break;}
idx = ((ix*ny+iy)*nz+iz)*nucell + iu;
tag2surf[itag] = idx;
surf2tag[idx] = itag;
}
fclose(fp);
if (tag2surf.size() != surf2tag.size() ||
if (tag2surf.size() != surf2tag.size() ||
tag2surf.size() != static_cast<std::size_t>(ngroup) )
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;
tagint *tag = atom->tag;
@ -616,7 +615,7 @@ void FixPhonon::readmap()
if (mask[i] & groupbit){
itag = tag[i];
idx = tag2surf[itag];
if (itag != surf2tag[idx])
if (itag != surf2tag[idx])
error->one(FLERR,"The mapping info read is incorrect!");
}
}
@ -631,11 +630,11 @@ void FixPhonon::postprocess( )
ifreq = 0;
int idim, jdim, ndim;
double inv_neval = 1.0 /double(neval);
double inv_neval = 1. /double(neval);
// to get <Rq.Rq*>
for (idq = 0; idq < mynq; ++idq)
for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim]*inv_neval;
for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim] * inv_neval;
// to get <R>
for (idx = 0; idx < mynpt; ++idx)
@ -645,15 +644,15 @@ void FixPhonon::postprocess( )
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_data[m++] = static_cast<FFT_SCALAR>(Rnow[idx][idim]);
fft_data[m++] = static_cast<FFT_SCALAR>(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]);
Rqnow[idq][idim] = std::complex<double>(static_cast<double>(fft_data[m]), static_cast<double>(fft_data[m+1]));
m += 2;
}
}
@ -662,20 +661,20 @@ void FixPhonon::postprocess( )
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]);
for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] -= Rqnow[idq][idim] * std::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 = inv_neval*inv_nTemp;
double NormFac = TempFac*double(ntotal);
double TempFac = inv_neval * inv_nTemp;
double NormFac = TempFac * double(ntotal);
for (idim = 0; idim < sysdim; ++idim){
kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac);
TempAve += TempSum[idim]*TempFac;
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;
@ -688,7 +687,7 @@ void FixPhonon::postprocess( )
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);
@ -707,7 +706,7 @@ void FixPhonon::postprocess( )
basevec[7] = hsum[3] * inv_neval / double(nz);
basevec[6] = hsum[4] * inv_neval / double(nz);
basevec[3] = hsum[5] * inv_neval / 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];
@ -822,7 +821,7 @@ void FixPhonon::GaussJordan(int n, std::complex<double> *Mat)
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;
@ -881,7 +880,7 @@ void FixPhonon::EnforceASR()
}
}
}
// symmetrize
for (int k = 0; k < nucell; ++k)
for (int kp = k; kp < nucell; ++kp){

View File

@ -31,9 +31,6 @@ FixStyle(phonon,FixPhonon)
#ifndef FIX_PHONON_H
#define FIX_PHONON_H
#include "lmptype.h"
#include <mpi.h>
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_FLOAT
@ -44,10 +41,11 @@ typedef double FFT_SCALAR;
#include <complex>
#include "fix.h"
#include <map>
#include <map.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "fft3d_wrap.h"
namespace LAMMPS_NS {
@ -75,7 +73,7 @@ class FixPhonon : public Fix {
int ngroup, nfind; // total number of atoms in group; total number of atoms on this proc
char *prefix, *logfile; // prefix of output file names
FILE *flog;
double *M_inv_sqrt;
class FFT3d *fft; // to do fft via the fft3d wraper
@ -84,7 +82,7 @@ class FixPhonon : public Fix {
int *fft_cnts, *fft_disp;
int fft_dim, fft_dim2;
FFT_SCALAR *fft_data;
tagint itag; // index variables
int idx, idq; // more index variables
std::map<tagint,int> tag2surf; // Mapping info

View File

@ -1163,7 +1163,6 @@ int ComputeChunkAtom::setup_xyz_bins()
// allocate and initialize arrays based on new bin count
double binlo[3],binhi[3];
double *prd;
if (scaleflag == REDUCED) {
binlo[0] = domain->boxlo_lamda[0];
binlo[1] = domain->boxlo_lamda[1];
@ -1171,7 +1170,6 @@ int ComputeChunkAtom::setup_xyz_bins()
binhi[0] = domain->boxhi_lamda[0];
binhi[1] = domain->boxhi_lamda[1];
binhi[2] = domain->boxhi_lamda[2];
prd = domain->prd_lamda;
} else {
binlo[0] = domain->boxlo[0];
binlo[1] = domain->boxlo[1];
@ -1179,7 +1177,6 @@ int ComputeChunkAtom::setup_xyz_bins()
binhi[0] = domain->boxhi[0];
binhi[1] = domain->boxhi[1];
binhi[2] = domain->boxhi[2];
prd = domain->prd;
}
if (minflag[0] == COORD) binlo[0] = minvalue[0];

View File

@ -305,7 +305,7 @@ void FixAdapt::init()
char *cptr;
int nsub = 0;
if (cptr = strchr(pstyle,':')) {
if ((cptr = strchr(pstyle,':'))) {
*cptr = '\0';
nsub = force->inumeric(FLERR,cptr+1);
}

View File

@ -197,6 +197,8 @@ typedef int bigint;
#define _noalias
#endif
#define ISFINITE(x) isfinite(x)
// settings to enable LAMMPS to build under Windows
#ifdef _WIN32