git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11907 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
@ -437,8 +437,7 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
|
|||||||
double qqrd2e = force->qqrd2e;
|
double qqrd2e = force->qqrd2e;
|
||||||
|
|
||||||
int i, j;
|
int i, j;
|
||||||
int order1 = ewald_order&(1<<1), order3 = ewald_order&(1<<3),
|
int order3 = ewald_order&(1<<3), order6 = ewald_order&(1<<6);
|
||||||
order6 = ewald_order&(1<<6);
|
|
||||||
int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
|
int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
|
||||||
double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
|
double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
|
||||||
double rsq, r2inv, force_coul, force_lj;
|
double rsq, r2inv, force_coul, force_lj;
|
||||||
|
|||||||
@ -80,13 +80,13 @@ PairLubricate::~PairLubricate()
|
|||||||
void PairLubricate::compute(int eflag, int vflag)
|
void PairLubricate::compute(int eflag, int vflag)
|
||||||
{
|
{
|
||||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||||
double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz;
|
double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
|
||||||
double rsq,r,h_sep,radi,tfmag;
|
double rsq,r,h_sep,radi;
|
||||||
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
|
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
|
||||||
double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
|
double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
|
||||||
double inertia,inv_inertia,vRS0;
|
double vRS0;
|
||||||
double vi[3],vj[3],wi[3],wj[3],xl[3];
|
double vi[3],vj[3],wi[3],wj[3],xl[3];
|
||||||
double a_sq,a_sh,a_pu,Fbmag,del,delmin,eta;
|
double a_sq,a_sh,a_pu;
|
||||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||||
double lamda[3],vstream[3];
|
double lamda[3],vstream[3];
|
||||||
|
|
||||||
@ -99,11 +99,8 @@ void PairLubricate::compute(int eflag, int vflag)
|
|||||||
double **v = atom->v;
|
double **v = atom->v;
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
double **omega = atom->omega;
|
double **omega = atom->omega;
|
||||||
double **angmom = atom->angmom;
|
|
||||||
double **torque = atom->torque;
|
double **torque = atom->torque;
|
||||||
double *radius = atom->radius;
|
double *radius = atom->radius;
|
||||||
double *mass = atom->mass;
|
|
||||||
double *rmass = atom->rmass;
|
|
||||||
int *type = atom->type;
|
int *type = atom->type;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
int newton_pair = force->newton_pair;
|
int newton_pair = force->newton_pair;
|
||||||
@ -546,7 +543,7 @@ void PairLubricate::init_style()
|
|||||||
if (comm->ghost_velocity == 0)
|
if (comm->ghost_velocity == 0)
|
||||||
error->all(FLERR,"Pair lubricate requires ghost atoms store velocity");
|
error->all(FLERR,"Pair lubricate requires ghost atoms store velocity");
|
||||||
|
|
||||||
int irequest = neighbor->request(this);
|
neighbor->request(this);
|
||||||
|
|
||||||
// require that atom radii are identical within each type
|
// require that atom radii are identical within each type
|
||||||
// require monodisperse system with same radii for all types
|
// require monodisperse system with same radii for all types
|
||||||
|
|||||||
@ -65,13 +65,13 @@ PairLubricatePoly::PairLubricatePoly(LAMMPS *lmp) : PairLubricate(lmp)
|
|||||||
void PairLubricatePoly::compute(int eflag, int vflag)
|
void PairLubricatePoly::compute(int eflag, int vflag)
|
||||||
{
|
{
|
||||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||||
double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz;
|
double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
|
||||||
double rsq,r,h_sep,h_sepj,beta0,beta1,betaj,betaj1,radi,radj,tfmag;
|
double rsq,r,h_sep,beta0,beta1,radi,radj;
|
||||||
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
|
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
|
||||||
double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
|
double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
|
||||||
double inertia,inv_inertia,vRS0;
|
double vRS0;
|
||||||
double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3];
|
double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3];
|
||||||
double a_sq,a_sh,a_pu,Fbmag,del,delmin,eta;
|
double a_sq,a_sh,a_pu;
|
||||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||||
double lamda[3],vstream[3];
|
double lamda[3],vstream[3];
|
||||||
|
|
||||||
@ -84,11 +84,8 @@ void PairLubricatePoly::compute(int eflag, int vflag)
|
|||||||
double **v = atom->v;
|
double **v = atom->v;
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
double **omega = atom->omega;
|
double **omega = atom->omega;
|
||||||
double **angmom = atom->angmom;
|
|
||||||
double **torque = atom->torque;
|
double **torque = atom->torque;
|
||||||
double *radius = atom->radius;
|
double *radius = atom->radius;
|
||||||
double *mass = atom->mass;
|
|
||||||
double *rmass = atom->rmass;
|
|
||||||
int *type = atom->type;
|
int *type = atom->type;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
int newton_pair = force->newton_pair;
|
int newton_pair = force->newton_pair;
|
||||||
@ -443,7 +440,6 @@ void PairLubricatePoly::init_style()
|
|||||||
// for pair hybrid, should limit test to types using the pair style
|
// for pair hybrid, should limit test to types using the pair style
|
||||||
|
|
||||||
double *radius = atom->radius;
|
double *radius = atom->radius;
|
||||||
int *type = atom->type;
|
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
for (int i = 0; i < nlocal; i++)
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
|||||||
@ -190,7 +190,7 @@ void PairCoulLongGPU::cpu_compute(int start, int inum, int eflag,
|
|||||||
int i,j,ii,jj,jnum,itable;
|
int i,j,ii,jj,jnum,itable;
|
||||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||||
double fraction,table;
|
double fraction,table;
|
||||||
double r,r2inv,r6inv,forcecoul,factor_coul;
|
double r,r2inv,forcecoul,factor_coul;
|
||||||
double grij,expm2,prefactor,t,erfc;
|
double grij,expm2,prefactor,t,erfc;
|
||||||
int *jlist;
|
int *jlist;
|
||||||
double rsq;
|
double rsq;
|
||||||
|
|||||||
@ -616,7 +616,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
|
|||||||
double radi,radj,radsum;
|
double radi,radj,radsum;
|
||||||
double r,rinv,rsqinv,delx,dely,delz;
|
double r,rinv,rsqinv,delx,dely,delz;
|
||||||
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
|
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
|
||||||
double mi,mj,meff,damp,ccel,polyhertz;
|
double mi,mj,meff,damp,ccel;
|
||||||
double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
|
double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
|
||||||
double fs1,fs2,fs3,fs,fn;
|
double fs1,fs2,fs3,fs,fn;
|
||||||
|
|
||||||
|
|||||||
@ -985,7 +985,6 @@ void MSM::set_grid_global()
|
|||||||
|
|
||||||
int flag = 0;
|
int flag = 0;
|
||||||
int xlevels,ylevels,zlevels;
|
int xlevels,ylevels,zlevels;
|
||||||
double k,r;
|
|
||||||
|
|
||||||
while (!factorable(nx_max,flag,xlevels)) {
|
while (!factorable(nx_max,flag,xlevels)) {
|
||||||
double k = log(nx_max)/log(2.0);
|
double k = log(nx_max)/log(2.0);
|
||||||
@ -1162,7 +1161,7 @@ void MSM::set_grid_local()
|
|||||||
|
|
||||||
// lengths of box and processor sub-domains
|
// lengths of box and processor sub-domains
|
||||||
|
|
||||||
double *prd,*sublo,*subhi,*boxhi;
|
double *prd,*sublo,*subhi;
|
||||||
|
|
||||||
if (!triclinic) {
|
if (!triclinic) {
|
||||||
prd = domain->prd;
|
prd = domain->prd;
|
||||||
@ -1546,7 +1545,7 @@ void MSM::direct(int n)
|
|||||||
|
|
||||||
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
||||||
int ii,jj,kk;
|
int ii,jj,kk;
|
||||||
int imin,imax,jmin,jmax,kmin,kmax;
|
int imin,imax,jmin,jmax,kmax;
|
||||||
double qtmp,qtmp2,gtmp;
|
double qtmp,qtmp2,gtmp;
|
||||||
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
|
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
|
||||||
double **qk,**ek;
|
double **qk,**ek;
|
||||||
@ -1560,10 +1559,8 @@ void MSM::direct(int n)
|
|||||||
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
||||||
|
|
||||||
if (domain->zperiodic) {
|
if (domain->zperiodic) {
|
||||||
kmin = nzlo_direct;
|
|
||||||
kmax = nzhi_direct;
|
kmax = nzhi_direct;
|
||||||
} else {
|
} else {
|
||||||
kmin = MAX(nzlo_direct,alpha[n] - icz);
|
|
||||||
kmax = MIN(nzhi_direct,betaz[n] - icz);
|
kmax = MIN(nzhi_direct,betaz[n] - icz);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -1759,7 +1756,7 @@ void MSM::direct_peratom(int n)
|
|||||||
|
|
||||||
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
||||||
int ii,jj,kk;
|
int ii,jj,kk;
|
||||||
int imin,imax,jmin,jmax,kmin,kmax;
|
int imin,imax,jmin,jmax,kmax;
|
||||||
double qtmp;
|
double qtmp;
|
||||||
|
|
||||||
int nx = nxhi_direct - nxlo_direct + 1;
|
int nx = nxhi_direct - nxlo_direct + 1;
|
||||||
@ -1770,10 +1767,8 @@ void MSM::direct_peratom(int n)
|
|||||||
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
||||||
|
|
||||||
if (domain->zperiodic) {
|
if (domain->zperiodic) {
|
||||||
kmin = nzlo_direct;
|
|
||||||
kmax = nzhi_direct;
|
kmax = nzhi_direct;
|
||||||
} else {
|
} else {
|
||||||
kmin = MAX(nzlo_direct,alpha[n] - icz);
|
|
||||||
kmax = MIN(nzhi_direct,betaz[n] - icz);
|
kmax = MIN(nzhi_direct,betaz[n] - icz);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -1901,7 +1896,7 @@ void MSM::direct_top(int n)
|
|||||||
|
|
||||||
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
||||||
int ii,jj,kk;
|
int ii,jj,kk;
|
||||||
int imin,imax,jmin,jmax,kmin,kmax;
|
int imin,imax,jmin,jmax,kmax;
|
||||||
double qtmp,qtmp2,gtmp;
|
double qtmp,qtmp2,gtmp;
|
||||||
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
|
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
|
||||||
double **qk,**ek;
|
double **qk,**ek;
|
||||||
@ -1919,10 +1914,8 @@ void MSM::direct_top(int n)
|
|||||||
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
||||||
|
|
||||||
if (domain->zperiodic) {
|
if (domain->zperiodic) {
|
||||||
kmin = 0;
|
|
||||||
kmax = nz_msm[n]-1;
|
kmax = nz_msm[n]-1;
|
||||||
} else {
|
} else {
|
||||||
kmin = alpha[n] - icz;
|
|
||||||
kmax = betaz[n] - icz;
|
kmax = betaz[n] - icz;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -2123,7 +2116,7 @@ void MSM::direct_peratom_top(int n)
|
|||||||
|
|
||||||
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
|
||||||
int ii,jj,kk;
|
int ii,jj,kk;
|
||||||
int imin,imax,jmin,jmax,kmin,kmax;
|
int imin,imax,jmin,jmax,kmax;
|
||||||
double qtmp;
|
double qtmp;
|
||||||
|
|
||||||
int nx_top = betax[n] - alpha[n];
|
int nx_top = betax[n] - alpha[n];
|
||||||
@ -2138,10 +2131,8 @@ void MSM::direct_peratom_top(int n)
|
|||||||
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
|
||||||
|
|
||||||
if (domain->zperiodic) {
|
if (domain->zperiodic) {
|
||||||
kmin = 0;
|
|
||||||
kmax = nz_msm[n]-1;
|
kmax = nz_msm[n]-1;
|
||||||
} else {
|
} else {
|
||||||
kmin = alpha[n] - icz;
|
|
||||||
kmax = betaz[n] - icz;
|
kmax = betaz[n] - icz;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
@ -404,7 +404,7 @@ void PPPMDisp::init()
|
|||||||
|
|
||||||
// print stats
|
// print stats
|
||||||
|
|
||||||
int ngrid_max,nfft_both_max,nbuf_max;
|
int ngrid_max,nfft_both_max;
|
||||||
MPI_Allreduce(&ngrid,&ngrid_max,1,MPI_INT,MPI_MAX,world);
|
MPI_Allreduce(&ngrid,&ngrid_max,1,MPI_INT,MPI_MAX,world);
|
||||||
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
|
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
|
||||||
|
|
||||||
@ -502,7 +502,7 @@ void PPPMDisp::init()
|
|||||||
|
|
||||||
// print stats
|
// print stats
|
||||||
|
|
||||||
int ngrid_max,nfft_both_max,nbuf_max;
|
int ngrid_max,nfft_both_max;
|
||||||
MPI_Allreduce(&ngrid_6,&ngrid_max,1,MPI_INT,MPI_MAX,world);
|
MPI_Allreduce(&ngrid_6,&ngrid_max,1,MPI_INT,MPI_MAX,world);
|
||||||
MPI_Allreduce(&nfft_both_6,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
|
MPI_Allreduce(&nfft_both_6,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
|
||||||
|
|
||||||
@ -2888,7 +2888,6 @@ double PPPMDisp::final_accuracy()
|
|||||||
|
|
||||||
void PPPMDisp::final_accuracy_6(double& acc, double& acc_real, double& acc_kspace)
|
void PPPMDisp::final_accuracy_6(double& acc, double& acc_real, double& acc_kspace)
|
||||||
{
|
{
|
||||||
double df_rspace, df_kspace;
|
|
||||||
double xprd = domain->xprd;
|
double xprd = domain->xprd;
|
||||||
double yprd = domain->yprd;
|
double yprd = domain->yprd;
|
||||||
double zprd = domain->zprd;
|
double zprd = domain->zprd;
|
||||||
@ -3050,12 +3049,10 @@ double PPPMDisp::compute_qopt_ad()
|
|||||||
double argx,argy,argz,wx,wy,wz,sx,sy,sz,qx,qy,qz;
|
double argx,argy,argz,wx,wy,wz,sx,sy,sz,qx,qy,qz;
|
||||||
double u2, sqk;
|
double u2, sqk;
|
||||||
double sum1,sum2,sum3,sum4,dot2;
|
double sum1,sum2,sum3,sum4,dot2;
|
||||||
double numerator;
|
|
||||||
|
|
||||||
int nbx = 2;
|
int nbx = 2;
|
||||||
int nby = 2;
|
int nby = 2;
|
||||||
int nbz = 2;
|
int nbz = 2;
|
||||||
double form = 1.0;
|
|
||||||
|
|
||||||
for (m = nzlo_fft; m <= nzhi_fft; m++) {
|
for (m = nzlo_fft; m <= nzhi_fft; m++) {
|
||||||
mper = m - nz_pppm*(2*m/nz_pppm);
|
mper = m - nz_pppm*(2*m/nz_pppm);
|
||||||
@ -3070,7 +3067,6 @@ double PPPMDisp::compute_qopt_ad()
|
|||||||
pow(unitkz*mper,2.0);
|
pow(unitkz*mper,2.0);
|
||||||
|
|
||||||
if (sqk != 0.0) {
|
if (sqk != 0.0) {
|
||||||
numerator = form*12.5663706;
|
|
||||||
|
|
||||||
sum1 = 0.0;
|
sum1 = 0.0;
|
||||||
sum2 = 0.0;
|
sum2 = 0.0;
|
||||||
@ -3120,7 +3116,7 @@ double PPPMDisp::compute_qopt_ad()
|
|||||||
double PPPMDisp::compute_qopt_6_ik()
|
double PPPMDisp::compute_qopt_6_ik()
|
||||||
{
|
{
|
||||||
double qopt = 0.0;
|
double qopt = 0.0;
|
||||||
int k,l,m,n;
|
int k,l,m;
|
||||||
double *prd;
|
double *prd;
|
||||||
|
|
||||||
if (triclinic == 0) prd = domain->prd;
|
if (triclinic == 0) prd = domain->prd;
|
||||||
@ -3148,7 +3144,6 @@ double PPPMDisp::compute_qopt_6_ik()
|
|||||||
int nby = 2;
|
int nby = 2;
|
||||||
int nbz = 2;
|
int nbz = 2;
|
||||||
|
|
||||||
n = 0;
|
|
||||||
for (m = nzlo_fft_6; m <= nzhi_fft_6; m++) {
|
for (m = nzlo_fft_6; m <= nzhi_fft_6; m++) {
|
||||||
mper = m - nz_pppm_6*(2*m/nz_pppm_6);
|
mper = m - nz_pppm_6*(2*m/nz_pppm_6);
|
||||||
|
|
||||||
@ -3343,9 +3338,6 @@ void PPPMDisp::calc_csum()
|
|||||||
//the following variables are needed to distinguish between arithmetic
|
//the following variables are needed to distinguish between arithmetic
|
||||||
// and geometric mixing
|
// and geometric mixing
|
||||||
|
|
||||||
double mix1; // scales 20/16 to 4
|
|
||||||
int mix2; // shifts the value to the sigma^3 value
|
|
||||||
int mix3; // shifts the value to the right atom type
|
|
||||||
if (function[1]) {
|
if (function[1]) {
|
||||||
for (i = 1; i <= ntypes; i++)
|
for (i = 1; i <= ntypes; i++)
|
||||||
cii[i] = B[i]*B[i];
|
cii[i] = B[i]*B[i];
|
||||||
|
|||||||
@ -263,7 +263,7 @@ void PPPMDispTIP4P::fieldforce_c_ik()
|
|||||||
void PPPMDispTIP4P::fieldforce_c_ad()
|
void PPPMDispTIP4P::fieldforce_c_ad()
|
||||||
{
|
{
|
||||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||||
FFT_SCALAR dx,dy,dz,x0,y0,z0,dx0,dy0,dz0;
|
FFT_SCALAR dx,dy,dz;
|
||||||
FFT_SCALAR ekx,eky,ekz;
|
FFT_SCALAR ekx,eky,ekz;
|
||||||
double *xi;
|
double *xi;
|
||||||
int iH1,iH2;
|
int iH1,iH2;
|
||||||
@ -399,7 +399,6 @@ void PPPMDispTIP4P::fieldforce_c_peratom()
|
|||||||
// ek = 3 components of E-field on particle
|
// ek = 3 components of E-field on particle
|
||||||
double *q = atom->q;
|
double *q = atom->q;
|
||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
double **f = atom->f;
|
|
||||||
|
|
||||||
int *type = atom->type;
|
int *type = atom->type;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|||||||
@ -278,7 +278,6 @@ void FixGCMC::init()
|
|||||||
|
|
||||||
if (molflag == 0 && atom->molecule_flag) {
|
if (molflag == 0 && atom->molecule_flag) {
|
||||||
tagint *molecule = atom->molecule;
|
tagint *molecule = atom->molecule;
|
||||||
int *mask = atom->mask;
|
|
||||||
int flag = 0;
|
int flag = 0;
|
||||||
for (int i = 0; i < atom->nlocal; i++)
|
for (int i = 0; i < atom->nlocal; i++)
|
||||||
if (type[i] == ngcmc_type)
|
if (type[i] == ngcmc_type)
|
||||||
@ -326,8 +325,8 @@ void FixGCMC::init()
|
|||||||
if (molflag) {
|
if (molflag) {
|
||||||
char **group_arg = new char*[3];
|
char **group_arg = new char*[3];
|
||||||
// create unique group name for atoms to be rotated
|
// create unique group name for atoms to be rotated
|
||||||
int len = strlen(id) + 24;
|
int len = strlen(id) + 30;
|
||||||
group_arg[0] = new char[60];
|
group_arg[0] = new char[len];
|
||||||
sprintf(group_arg[0],"FixGCMC:rotation_gas_atoms:%s",id);
|
sprintf(group_arg[0],"FixGCMC:rotation_gas_atoms:%s",id);
|
||||||
group_arg[1] = (char *) "molecule";
|
group_arg[1] = (char *) "molecule";
|
||||||
char digits[12];
|
char digits[12];
|
||||||
@ -778,10 +777,6 @@ void FixGCMC::attempt_molecule_insertion()
|
|||||||
{
|
{
|
||||||
ninsertion_attempts += 1.0;
|
ninsertion_attempts += 1.0;
|
||||||
|
|
||||||
double xprd = domain->xprd;
|
|
||||||
double yprd = domain->yprd;
|
|
||||||
double zprd = domain->zprd;
|
|
||||||
|
|
||||||
double com_coord[3];
|
double com_coord[3];
|
||||||
if (regionflag) {
|
if (regionflag) {
|
||||||
int region_attempt = 0;
|
int region_attempt = 0;
|
||||||
@ -857,7 +852,6 @@ void FixGCMC::attempt_molecule_insertion()
|
|||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
double **v = atom->v;
|
double **v = atom->v;
|
||||||
imageint *image = atom->image;
|
imageint *image = atom->image;
|
||||||
tagint *molecule = atom->molecule;
|
|
||||||
tagint *tag = atom->tag;
|
tagint *tag = atom->tag;
|
||||||
for (int i = 0; i < natoms_per_molecule; i++) {
|
for (int i = 0; i < natoms_per_molecule; i++) {
|
||||||
k += atom->avec->unpack_exchange(&model_atom_buf[k]);
|
k += atom->avec->unpack_exchange(&model_atom_buf[k]);
|
||||||
|
|||||||
@ -669,12 +669,10 @@ void FixRigid::init()
|
|||||||
void FixRigid::setup(int vflag)
|
void FixRigid::setup(int vflag)
|
||||||
{
|
{
|
||||||
int i,n,ibody;
|
int i,n,ibody;
|
||||||
double massone,radone;
|
|
||||||
|
|
||||||
// fcm = force on center-of-mass of each rigid body
|
// fcm = force on center-of-mass of each rigid body
|
||||||
|
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
int *type = atom->type;
|
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
for (ibody = 0; ibody < nbody; ibody++)
|
for (ibody = 0; ibody < nbody; ibody++)
|
||||||
@ -2032,7 +2030,7 @@ void FixRigid::setup_bodies_static()
|
|||||||
|
|
||||||
void FixRigid::setup_bodies_dynamic()
|
void FixRigid::setup_bodies_dynamic()
|
||||||
{
|
{
|
||||||
int i,n,ibody;
|
int i,ibody;
|
||||||
double massone,radone;
|
double massone,radone;
|
||||||
|
|
||||||
// vcm = velocity of center-of-mass of each rigid body
|
// vcm = velocity of center-of-mass of each rigid body
|
||||||
@ -2131,7 +2129,7 @@ void FixRigid::setup_bodies_dynamic()
|
|||||||
|
|
||||||
void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
|
void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
|
||||||
{
|
{
|
||||||
int i,j,m,nchunk,id,eofflag;
|
int j,nchunk,id,eofflag;
|
||||||
int nlines;
|
int nlines;
|
||||||
FILE *fp;
|
FILE *fp;
|
||||||
char *eof,*start,*next,*buf;
|
char *eof,*start,*next,*buf;
|
||||||
|
|||||||
@ -306,11 +306,9 @@ void FixRigidNHSmall::setup(int vflag)
|
|||||||
// total translational and rotational degrees of freedom
|
// total translational and rotational degrees of freedom
|
||||||
|
|
||||||
int k,ibody;
|
int k,ibody;
|
||||||
double *inertia;
|
|
||||||
|
|
||||||
nf_t = nf_r = dimension * nlocal_body;
|
nf_t = nf_r = dimension * nlocal_body;
|
||||||
for (ibody = 0; ibody < nlocal_body; ibody++) {
|
for (ibody = 0; ibody < nlocal_body; ibody++) {
|
||||||
inertia = body[ibody].inertia;
|
|
||||||
for (k = 0; k < domain->dimension; k++)
|
for (k = 0; k < domain->dimension; k++)
|
||||||
if (fabs(body[ibody].inertia[k]) < EPSILON) nf_r--;
|
if (fabs(body[ibody].inertia[k]) < EPSILON) nf_r--;
|
||||||
}
|
}
|
||||||
@ -625,7 +623,7 @@ void FixRigidNHSmall::final_integrate()
|
|||||||
{
|
{
|
||||||
int i,ibody;
|
int i,ibody;
|
||||||
double tmp,scale_t[3],scale_r;
|
double tmp,scale_t[3],scale_r;
|
||||||
double dtfm,xy,xz,yz;
|
double dtfm;
|
||||||
double mbody[3],tbody[3],fquat[4];
|
double mbody[3],tbody[3],fquat[4];
|
||||||
double dtf2 = dtf * 2.0;
|
double dtf2 = dtf * 2.0;
|
||||||
|
|
||||||
@ -1001,11 +999,11 @@ void FixRigidNHSmall::nhc_press_integrate()
|
|||||||
|
|
||||||
double FixRigidNHSmall::compute_scalar()
|
double FixRigidNHSmall::compute_scalar()
|
||||||
{
|
{
|
||||||
int i,k,ibody;
|
int i,k;
|
||||||
double kt = boltz * t_target;
|
double kt = boltz * t_target;
|
||||||
double energy,ke_t,ke_q,tmp,Pkq[4];
|
double energy,ke_t,ke_q,tmp,Pkq[4];
|
||||||
|
|
||||||
double *vcm,*inertia,*quat;
|
double *vcm,*quat;
|
||||||
|
|
||||||
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
|
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
|
||||||
|
|
||||||
|
|||||||
@ -299,7 +299,6 @@ void FixMSST::init()
|
|||||||
double mass = 0.0;
|
double mass = 0.0;
|
||||||
for (int i = 0; i < atom->nlocal; i++) mass += atom->mass[atom->type[i]];
|
for (int i = 0; i < atom->nlocal; i++) mass += atom->mass[atom->type[i]];
|
||||||
MPI_Allreduce(&mass,&total_mass,1,MPI_DOUBLE,MPI_SUM,world);
|
MPI_Allreduce(&mass,&total_mass,1,MPI_DOUBLE,MPI_SUM,world);
|
||||||
total_mass = total_mass;
|
|
||||||
|
|
||||||
if (force->kspace) kspace_flag = 1;
|
if (force->kspace) kspace_flag = 1;
|
||||||
else kspace_flag = 0;
|
else kspace_flag = 0;
|
||||||
|
|||||||
@ -68,7 +68,6 @@ FixLbPC::FixLbPC(LAMMPS *lmp, int narg, char **arg) :
|
|||||||
|
|
||||||
int *mask = atom->mask;
|
int *mask = atom->mask;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
int print_warning = 0;
|
|
||||||
for(int j=0; j<nlocal; j++){
|
for(int j=0; j<nlocal; j++){
|
||||||
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
|
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
|
||||||
error->one(FLERR,"can only use the lb/pc fix for an atom if also using the lb/fluid fix for that atom");
|
error->one(FLERR,"can only use the lb/pc fix for an atom if also using the lb/fluid fix for that atom");
|
||||||
|
|||||||
@ -455,7 +455,7 @@ int FixLbRigidPCSphere::setmask()
|
|||||||
|
|
||||||
void FixLbRigidPCSphere::init()
|
void FixLbRigidPCSphere::init()
|
||||||
{
|
{
|
||||||
int i,itype,ibody;
|
int i,ibody;
|
||||||
|
|
||||||
// warn if more than one rigid fix
|
// warn if more than one rigid fix
|
||||||
|
|
||||||
@ -649,7 +649,7 @@ void FixLbRigidPCSphere::init()
|
|||||||
void FixLbRigidPCSphere::setup(int vflag)
|
void FixLbRigidPCSphere::setup(int vflag)
|
||||||
{
|
{
|
||||||
int i,n,ibody;
|
int i,n,ibody;
|
||||||
double massone,radone;
|
double massone;
|
||||||
|
|
||||||
// vcm = velocity of center-of-mass of each rigid body
|
// vcm = velocity of center-of-mass of each rigid body
|
||||||
// fcm = force on center-of-mass of each rigid body
|
// fcm = force on center-of-mass of each rigid body
|
||||||
@ -663,7 +663,6 @@ void FixLbRigidPCSphere::setup(int vflag)
|
|||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
imageint *image = atom->image;
|
imageint *image = atom->image;
|
||||||
int *periodicity = domain->periodicity;
|
|
||||||
|
|
||||||
double unwrap[3];
|
double unwrap[3];
|
||||||
double dx,dy,dz;
|
double dx,dy,dz;
|
||||||
@ -761,12 +760,11 @@ void FixLbRigidPCSphere::initial_integrate(int vflag)
|
|||||||
{
|
{
|
||||||
double dtfm;
|
double dtfm;
|
||||||
|
|
||||||
int i,n,ibody;
|
int i,ibody;
|
||||||
|
|
||||||
double massone,radone;
|
double massone;
|
||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
double **v = atom->v;
|
double **v = atom->v;
|
||||||
double **f = atom->f;
|
|
||||||
double *rmass = atom->rmass;
|
double *rmass = atom->rmass;
|
||||||
double *mass = atom->mass;
|
double *mass = atom->mass;
|
||||||
int *type = atom->type;
|
int *type = atom->type;
|
||||||
@ -953,14 +951,12 @@ void FixLbRigidPCSphere::initial_integrate(int vflag)
|
|||||||
void FixLbRigidPCSphere::final_integrate()
|
void FixLbRigidPCSphere::final_integrate()
|
||||||
{
|
{
|
||||||
int i,ibody;
|
int i,ibody;
|
||||||
double xy,xz,yz;
|
|
||||||
|
|
||||||
// sum over atoms to get force and torque on rigid body
|
// sum over atoms to get force and torque on rigid body
|
||||||
double massone;
|
double massone;
|
||||||
imageint *image = atom->image;
|
imageint *image = atom->image;
|
||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
double **v = atom->v;
|
|
||||||
double *rmass = atom->rmass;
|
double *rmass = atom->rmass;
|
||||||
double *mass = atom->mass;
|
double *mass = atom->mass;
|
||||||
int *type = atom->type;
|
int *type = atom->type;
|
||||||
@ -1565,7 +1561,7 @@ double FixLbRigidPCSphere::compute_array(int i, int j)
|
|||||||
int *mask = atom->mask;
|
int *mask = atom->mask;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
int i,j,k,p,m;
|
int i,k;
|
||||||
int ix,iy,iz;
|
int ix,iy,iz;
|
||||||
int ixp,iyp,izp;
|
int ixp,iyp,izp;
|
||||||
double dx1,dy1,dz1;
|
double dx1,dy1,dz1;
|
||||||
|
|||||||
@ -299,11 +299,11 @@ void PairCDEAM::compute(int eflag, int vflag)
|
|||||||
// Calculate derivative of h(x_i) polynomial function.
|
// Calculate derivative of h(x_i) polynomial function.
|
||||||
h_prime_i = evalHprime(x_i);
|
h_prime_i = evalHprime(x_i);
|
||||||
D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
|
D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
|
||||||
}
|
} else if(cdeamVersion == 2) {
|
||||||
else if(cdeamVersion == 2) {
|
|
||||||
D_i = D_values[i];
|
D_i = D_values[i];
|
||||||
|
} else {
|
||||||
|
ASSERT(false);
|
||||||
}
|
}
|
||||||
else ASSERT(false);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
for(jj = 0; jj < jnum; jj++) {
|
for(jj = 0; jj < jnum; jj++) {
|
||||||
@ -347,12 +347,11 @@ void PairCDEAM::compute(int eflag, int vflag)
|
|||||||
// Calculate derivative of h(x_j) polynomial function.
|
// Calculate derivative of h(x_j) polynomial function.
|
||||||
double h_prime_j = evalHprime(x_j);
|
double h_prime_j = evalHprime(x_j);
|
||||||
D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
|
D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
|
||||||
}
|
} else if(cdeamVersion == 2) {
|
||||||
else if(cdeamVersion == 2) {
|
|
||||||
D_j = D_values[j];
|
D_j = D_values[j];
|
||||||
|
} else {
|
||||||
|
ASSERT(false);
|
||||||
}
|
}
|
||||||
else ASSERT(false);
|
|
||||||
|
|
||||||
double t2 = -rhoB[j];
|
double t2 = -rhoB[j];
|
||||||
if(itype == speciesB) t2 += rho[j];
|
if(itype == speciesB) t2 += rho[j];
|
||||||
fpair += D_j * rhoip * t2;
|
fpair += D_j * rhoip * t2;
|
||||||
@ -372,8 +371,7 @@ void PairCDEAM::compute(int eflag, int vflag)
|
|||||||
if(itype == jtype || x_i == -1.0 || x_j == -1.0) {
|
if(itype == jtype || x_i == -1.0 || x_j == -1.0) {
|
||||||
// Case of no concentration dependence.
|
// Case of no concentration dependence.
|
||||||
fpair += phip;
|
fpair += phip;
|
||||||
}
|
} else {
|
||||||
else {
|
|
||||||
// We have a concentration dependence for the i-j interaction.
|
// We have a concentration dependence for the i-j interaction.
|
||||||
double h=0.0;
|
double h=0.0;
|
||||||
if(cdeamVersion == 1) {
|
if(cdeamVersion == 1) {
|
||||||
@ -382,14 +380,14 @@ void PairCDEAM::compute(int eflag, int vflag)
|
|||||||
// Calculate h(x_j) polynomial function.
|
// Calculate h(x_j) polynomial function.
|
||||||
double h_j = evalH(x_j);
|
double h_j = evalH(x_j);
|
||||||
h = 0.5 * (h_i + h_j);
|
h = 0.5 * (h_i + h_j);
|
||||||
}
|
} else if(cdeamVersion == 2) {
|
||||||
else if(cdeamVersion == 2) {
|
|
||||||
// Average concentration.
|
// Average concentration.
|
||||||
double x_ij = 0.5 * (x_i + x_j);
|
double x_ij = 0.5 * (x_i + x_j);
|
||||||
// Calculate h(x_ij) polynomial function.
|
// Calculate h(x_ij) polynomial function.
|
||||||
h = evalH(x_ij);
|
h = evalH(x_ij);
|
||||||
|
} else {
|
||||||
|
ASSERT(false);
|
||||||
}
|
}
|
||||||
else ASSERT(false);
|
|
||||||
fpair += h * phip;
|
fpair += h * phip;
|
||||||
phi *= h;
|
phi *= h;
|
||||||
}
|
}
|
||||||
@ -556,8 +554,9 @@ void PairCDEAM::unpack_comm(int n, int first, double *buf)
|
|||||||
rho[i] = buf[m++];
|
rho[i] = buf[m++];
|
||||||
rhoB[i] = buf[m++];
|
rhoB[i] = buf[m++];
|
||||||
}
|
}
|
||||||
|
} else {
|
||||||
|
ASSERT(false);
|
||||||
}
|
}
|
||||||
else ASSERT(false);
|
|
||||||
}
|
}
|
||||||
else if(communicationStage == 4) {
|
else if(communicationStage == 4) {
|
||||||
for(i = first; i < last; i++) {
|
for(i = first; i < last; i++) {
|
||||||
@ -616,15 +615,15 @@ void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf)
|
|||||||
rhoB[j] += buf[m++];
|
rhoB[j] += buf[m++];
|
||||||
D_values[j] += buf[m++];
|
D_values[j] += buf[m++];
|
||||||
}
|
}
|
||||||
}
|
} else if(cdeamVersion == 2) {
|
||||||
else if(cdeamVersion == 2) {
|
|
||||||
for(i = 0; i < n; i++) {
|
for(i = 0; i < n; i++) {
|
||||||
j = list[i];
|
j = list[i];
|
||||||
rho[j] += buf[m++];
|
rho[j] += buf[m++];
|
||||||
rhoB[j] += buf[m++];
|
rhoB[j] += buf[m++];
|
||||||
}
|
}
|
||||||
|
} else {
|
||||||
|
ASSERT(false);
|
||||||
}
|
}
|
||||||
else ASSERT(false);
|
|
||||||
}
|
}
|
||||||
else if(communicationStage == 3) {
|
else if(communicationStage == 3) {
|
||||||
for(i = 0; i < n; i++) {
|
for(i = 0; i < n; i++) {
|
||||||
|
|||||||
@ -91,7 +91,7 @@ PairEDIP::~PairEDIP()
|
|||||||
void PairEDIP::compute(int eflag, int vflag)
|
void PairEDIP::compute(int eflag, int vflag)
|
||||||
{
|
{
|
||||||
int i,j,k,ii,inum,jnum;
|
int i,j,k,ii,inum,jnum;
|
||||||
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
|
int itype,jtype,ktype,ijparam,ikparam;
|
||||||
double xtmp,ytmp,ztmp,evdwl;
|
double xtmp,ytmp,ztmp,evdwl;
|
||||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||||
register int preForceCoord_counter;
|
register int preForceCoord_counter;
|
||||||
@ -127,7 +127,6 @@ void PairEDIP::compute(int eflag, int vflag)
|
|||||||
double exp3B_ik;
|
double exp3B_ik;
|
||||||
double exp3BDerived_ik;
|
double exp3BDerived_ik;
|
||||||
double qFunction;
|
double qFunction;
|
||||||
double qFunctionDerived;
|
|
||||||
double tauFunction;
|
double tauFunction;
|
||||||
double tauFunctionDerived;
|
double tauFunctionDerived;
|
||||||
double expMinusBetaZeta_iZeta_i;
|
double expMinusBetaZeta_iZeta_i;
|
||||||
@ -290,8 +289,6 @@ void PairEDIP::compute(int eflag, int vflag)
|
|||||||
tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) *
|
tauFunctionDerived = interpolY1 + (interpolY2 - interpolY1) *
|
||||||
(interpolTMP-interpolIDX);
|
(interpolTMP-interpolIDX);
|
||||||
|
|
||||||
qFunctionDerived = -mu * qFunction;
|
|
||||||
|
|
||||||
forceModCoord_factor = 2.0 * beta * zeta_i * expMinusBetaZeta_iZeta_i;
|
forceModCoord_factor = 2.0 * beta * zeta_i * expMinusBetaZeta_iZeta_i;
|
||||||
|
|
||||||
forceModCoord = 0.0;
|
forceModCoord = 0.0;
|
||||||
@ -365,7 +362,6 @@ void PairEDIP::compute(int eflag, int vflag)
|
|||||||
k &= NEIGHMASK;
|
k &= NEIGHMASK;
|
||||||
ktype = map[type[k]];
|
ktype = map[type[k]];
|
||||||
ikparam = elem2param[itype][ktype][ktype];
|
ikparam = elem2param[itype][ktype][ktype];
|
||||||
ijkparam = elem2param[itype][jtype][ktype];
|
|
||||||
|
|
||||||
dr_ik[0] = x[k][0] - xtmp;
|
dr_ik[0] = x[k][0] - xtmp;
|
||||||
dr_ik[1] = x[k][1] - ytmp;
|
dr_ik[1] = x[k][1] - ytmp;
|
||||||
@ -639,7 +635,6 @@ void PairEDIP::initGrids(void)
|
|||||||
int numGridPointsNotOneCutoffFunction;
|
int numGridPointsNotOneCutoffFunction;
|
||||||
int numGridPointsCutoffFunction;
|
int numGridPointsCutoffFunction;
|
||||||
int numGridPointsR;
|
int numGridPointsR;
|
||||||
int numGridPointsRTotal;
|
|
||||||
int numGridPointsQFunctionGrid;
|
int numGridPointsQFunctionGrid;
|
||||||
int numGridPointsExpMinusBetaZeta_iZeta_i;
|
int numGridPointsExpMinusBetaZeta_iZeta_i;
|
||||||
int numGridPointsTauFunctionGrid;
|
int numGridPointsTauFunctionGrid;
|
||||||
@ -737,7 +732,6 @@ void PairEDIP::initGrids(void)
|
|||||||
|
|
||||||
numGridPointsR = (int)
|
numGridPointsR = (int)
|
||||||
((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY);
|
((cutoffA + leftLimitToZero - GRIDSTART) * GRIDDENSITY);
|
||||||
numGridPointsRTotal = numGridPointsR + 2;
|
|
||||||
|
|
||||||
r = GRIDSTART;
|
r = GRIDSTART;
|
||||||
deltaArgumentR = 1.0 / GRIDDENSITY;
|
deltaArgumentR = 1.0 / GRIDDENSITY;
|
||||||
|
|||||||
@ -84,13 +84,13 @@ void PairNMCutCoulLongOMP::compute(int eflag, int vflag)
|
|||||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||||
void PairNMCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
void PairNMCutCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||||
{
|
{
|
||||||
int i,j,ii,jj,jnum,jtype,itable;
|
int j,ii,jj,jnum,jtype,itable;
|
||||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||||
double fraction,table;
|
double fraction,table;
|
||||||
double r,rsq,rinv,r2inv,factor_coul,factor_lj;
|
double r,rsq,r2inv,factor_coul,factor_lj;
|
||||||
double forcecoul,forcenm,rminv,rninv;
|
double forcecoul,forcenm,rminv,rninv;
|
||||||
double grij,expm2,prefactor,t,erfc;
|
double grij,expm2,prefactor,t,erfc;
|
||||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
int *ilist,*numneigh,**firstneigh;
|
||||||
|
|
||||||
evdwl = ecoul = 0.0;
|
evdwl = ecoul = 0.0;
|
||||||
|
|
||||||
|
|||||||
@ -775,7 +775,6 @@ void PPPMDispTIP4POMP::make_rho_a()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_c_ik()
|
void PPPMDispTIP4POMP::fieldforce_c_ik()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -806,7 +805,7 @@ void PPPMDispTIP4POMP::fieldforce_c_ik()
|
|||||||
FFT_SCALAR x0,y0,z0,ekx,eky,ekz;
|
FFT_SCALAR x0,y0,z0,ekx,eky,ekz;
|
||||||
int i,ifrom,ito,tid,iH1,iH2,l,m,n,mx,my,mz;
|
int i,ifrom,ito,tid,iH1,iH2,l,m,n,mx,my,mz;
|
||||||
|
|
||||||
loop_setup_thr(ifrom,ito,tid,nlocal,nthreads);
|
loop_setup_thr(ifrom,ito,tid,nlocal,comm->nthreads);
|
||||||
|
|
||||||
// get per thread data
|
// get per thread data
|
||||||
ThrData *thr = fix->get_thr(tid);
|
ThrData *thr = fix->get_thr(tid);
|
||||||
@ -879,7 +878,6 @@ void PPPMDispTIP4POMP::fieldforce_c_ik()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_c_ad()
|
void PPPMDispTIP4POMP::fieldforce_c_ad()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -916,7 +914,7 @@ void PPPMDispTIP4POMP::fieldforce_c_ad()
|
|||||||
FFT_SCALAR ekx,eky,ekz;
|
FFT_SCALAR ekx,eky,ekz;
|
||||||
int i,ifrom,ito,tid,iH1,iH2,l,m,n,mx,my,mz;
|
int i,ifrom,ito,tid,iH1,iH2,l,m,n,mx,my,mz;
|
||||||
|
|
||||||
loop_setup_thr(ifrom,ito,tid,nlocal,nthreads);
|
loop_setup_thr(ifrom,ito,tid,nlocal,comm->nthreads);
|
||||||
|
|
||||||
// get per thread data
|
// get per thread data
|
||||||
ThrData *thr = fix->get_thr(tid);
|
ThrData *thr = fix->get_thr(tid);
|
||||||
@ -1008,7 +1006,6 @@ void PPPMDispTIP4POMP::fieldforce_c_ad()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_g_ik()
|
void PPPMDispTIP4POMP::fieldforce_g_ik()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1031,7 +1028,7 @@ void PPPMDispTIP4POMP::fieldforce_g_ik()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
@ -1097,7 +1094,6 @@ void PPPMDispTIP4POMP::fieldforce_g_ik()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_g_ad()
|
void PPPMDispTIP4POMP::fieldforce_g_ad()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1133,7 +1129,7 @@ void PPPMDispTIP4POMP::fieldforce_g_ad()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
@ -1147,7 +1143,7 @@ void PPPMDispTIP4POMP::fieldforce_g_ad()
|
|||||||
FFT_SCALAR * const * const dr1d = static_cast<FFT_SCALAR **>(thr->get_drho1d_6());
|
FFT_SCALAR * const * const dr1d = static_cast<FFT_SCALAR **>(thr->get_drho1d_6());
|
||||||
|
|
||||||
int l,m,n,nx,ny,nz,mx,my,mz;
|
int l,m,n,nx,ny,nz,mx,my,mz;
|
||||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
FFT_SCALAR dx,dy,dz;
|
||||||
FFT_SCALAR ekx,eky,ekz;
|
FFT_SCALAR ekx,eky,ekz;
|
||||||
int type;
|
int type;
|
||||||
double lj;
|
double lj;
|
||||||
@ -1219,7 +1215,6 @@ void PPPMDispTIP4POMP::fieldforce_g_ad()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_g_peratom()
|
void PPPMDispTIP4POMP::fieldforce_g_peratom()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1241,7 +1236,7 @@ void PPPMDispTIP4POMP::fieldforce_g_peratom()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
@ -1252,7 +1247,7 @@ void PPPMDispTIP4POMP::fieldforce_g_peratom()
|
|||||||
ThrData *thr = fix->get_thr(tid);
|
ThrData *thr = fix->get_thr(tid);
|
||||||
FFT_SCALAR * const * const r1d = static_cast<FFT_SCALAR **>(thr->get_rho1d_6());
|
FFT_SCALAR * const * const r1d = static_cast<FFT_SCALAR **>(thr->get_rho1d_6());
|
||||||
|
|
||||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
int l,m,n,nx,ny,nz,mx,my,mz;
|
||||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
||||||
FFT_SCALAR u,v0,v1,v2,v3,v4,v5;
|
FFT_SCALAR u,v0,v1,v2,v3,v4,v5;
|
||||||
int type;
|
int type;
|
||||||
@ -1318,7 +1313,6 @@ void PPPMDispTIP4POMP::fieldforce_g_peratom()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_a_ik()
|
void PPPMDispTIP4POMP::fieldforce_a_ik()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1332,7 +1326,6 @@ void PPPMDispTIP4POMP::fieldforce_a_ik()
|
|||||||
// ek = 3 components of E-field on particle
|
// ek = 3 components of E-field on particle
|
||||||
|
|
||||||
const double * const * const x = atom->x;
|
const double * const * const x = atom->x;
|
||||||
const double qqrd2e = force->qqrd2e;
|
|
||||||
|
|
||||||
#if defined(_OPENMP)
|
#if defined(_OPENMP)
|
||||||
#pragma omp parallel default(none)
|
#pragma omp parallel default(none)
|
||||||
@ -1342,7 +1335,7 @@ void PPPMDispTIP4POMP::fieldforce_a_ik()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
@ -1440,7 +1433,6 @@ void PPPMDispTIP4POMP::fieldforce_a_ik()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_a_ad()
|
void PPPMDispTIP4POMP::fieldforce_a_ad()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1476,7 +1468,7 @@ void PPPMDispTIP4POMP::fieldforce_a_ad()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
@ -1629,7 +1621,6 @@ void PPPMDispTIP4POMP::fieldforce_a_ad()
|
|||||||
|
|
||||||
void PPPMDispTIP4POMP::fieldforce_a_peratom()
|
void PPPMDispTIP4POMP::fieldforce_a_peratom()
|
||||||
{
|
{
|
||||||
const int nthreads = comm->nthreads;
|
|
||||||
const int nlocal = atom->nlocal;
|
const int nlocal = atom->nlocal;
|
||||||
|
|
||||||
// no local atoms => nothing to do
|
// no local atoms => nothing to do
|
||||||
@ -1651,7 +1642,7 @@ void PPPMDispTIP4POMP::fieldforce_a_peratom()
|
|||||||
// each thread works on a fixed chunk of atoms.
|
// each thread works on a fixed chunk of atoms.
|
||||||
const int tid = omp_get_thread_num();
|
const int tid = omp_get_thread_num();
|
||||||
const int inum = nlocal;
|
const int inum = nlocal;
|
||||||
const int idelta = 1 + inum/nthreads;
|
const int idelta = 1 + inum/comm->nthreads;
|
||||||
const int ifrom = tid*idelta;
|
const int ifrom = tid*idelta;
|
||||||
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
|
||||||
#else
|
#else
|
||||||
|
|||||||
@ -170,8 +170,7 @@ void ComputeVoronoi::compute_peratom()
|
|||||||
// decide between occupation or per-frame tesselation modes
|
// decide between occupation or per-frame tesselation modes
|
||||||
if (occupation) {
|
if (occupation) {
|
||||||
// build cells only once
|
// build cells only once
|
||||||
int i, j,
|
int i, nall = nlocal + atom->nghost;
|
||||||
nall = nlocal + atom->nghost;
|
|
||||||
if (con_mono==NULL && con_poly==NULL) {
|
if (con_mono==NULL && con_poly==NULL) {
|
||||||
// generate the voronoi cell network for the initial structure
|
// generate the voronoi cell network for the initial structure
|
||||||
buildCells();
|
buildCells();
|
||||||
@ -205,7 +204,7 @@ void ComputeVoronoi::compute_peratom()
|
|||||||
|
|
||||||
void ComputeVoronoi::buildCells()
|
void ComputeVoronoi::buildCells()
|
||||||
{
|
{
|
||||||
int i, j;
|
int i;
|
||||||
const double e = 0.01;
|
const double e = 0.01;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
@ -219,7 +218,7 @@ void ComputeVoronoi::buildCells()
|
|||||||
}
|
}
|
||||||
|
|
||||||
double *sublo = domain->sublo, *sublo_lamda = domain->sublo_lamda, *boxlo = domain->boxlo;
|
double *sublo = domain->sublo, *sublo_lamda = domain->sublo_lamda, *boxlo = domain->boxlo;
|
||||||
double *subhi = domain->subhi, *subhi_lamda = domain->subhi_lamda, *boxhi = domain->boxhi;
|
double *subhi = domain->subhi, *subhi_lamda = domain->subhi_lamda;
|
||||||
double *cut = comm->cutghost;
|
double *cut = comm->cutghost;
|
||||||
double sublo_bound[3], subhi_bound[3], cut_bound[3];
|
double sublo_bound[3], subhi_bound[3], cut_bound[3];
|
||||||
double **x = atom->x;
|
double **x = atom->x;
|
||||||
@ -229,7 +228,7 @@ void ComputeVoronoi::buildCells()
|
|||||||
// triclinic box: embed parallelepiped into orthogonal voro++ domain
|
// triclinic box: embed parallelepiped into orthogonal voro++ domain
|
||||||
|
|
||||||
// cutghost is in lamda coordinates for triclinic boxes, use subxx_lamda
|
// cutghost is in lamda coordinates for triclinic boxes, use subxx_lamda
|
||||||
double *h = domain->h, cuttri[3];
|
double *h = domain->h;
|
||||||
sublo_bound[0] = h[0]*sublo_lamda[0] + h[5]*sublo_lamda[1] + h[4]*sublo_lamda[2] + boxlo[0];
|
sublo_bound[0] = h[0]*sublo_lamda[0] + h[5]*sublo_lamda[1] + h[4]*sublo_lamda[2] + boxlo[0];
|
||||||
sublo_bound[1] = h[1]*sublo_lamda[1] + h[3]*sublo_lamda[2] + boxlo[1];
|
sublo_bound[1] = h[1]*sublo_lamda[1] + h[3]*sublo_lamda[2] + boxlo[1];
|
||||||
sublo_bound[2] = h[2]*sublo_lamda[2] + boxlo[2];
|
sublo_bound[2] = h[2]*sublo_lamda[2] + boxlo[2];
|
||||||
|
|||||||
@ -72,6 +72,9 @@ void WriteDump::command(int narg, char **arg)
|
|||||||
if (strcmp(arg[1],"image") == 0)
|
if (strcmp(arg[1],"image") == 0)
|
||||||
((DumpImage *) dump)->multifile_override = 1;
|
((DumpImage *) dump)->multifile_override = 1;
|
||||||
|
|
||||||
|
if (strcmp(arg[1],"cfg") == 0)
|
||||||
|
((DumpCFG *) dump)->multifile_override = 1;
|
||||||
|
|
||||||
dump->init();
|
dump->init();
|
||||||
dump->write();
|
dump->write();
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user