jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

whitespace cleanup: remove trailing blanks

This commit is contained in:
Axel Kohlmeyer 2019-11-03 11:03:39 -05:00
parent 9e7ca428aa
commit b6b022b610
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
69 changed files with 452 additions and 452 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/CLASS2/pair_lj_class2_coul_long.cpp
View File

@ -664,19 +664,19 @@ void PairLJClass2CoulLong::init_style()
if (!atom->q_flag)
error->all(FLERR,
"Pair style lj/class2/coul/long requires atom attribute q");
// request regular or rRESPA neighbor list
int irequest;
int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
irequest = neighbor->request(this,instance_me);
if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1;
neighbor->requests[irequest]->respainner = 1;
@ -684,13 +684,13 @@ void PairLJClass2CoulLong::init_style()
if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
cut_coulsq = cut_coul * cut_coul;
// set rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = NULL;
else cut_respa = NULL;
// insure use of KSpace long-range solver, set g_ewald
@ -739,9 +739,9 @@ double PairLJClass2CoulLong::init_one(int i, int j)
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
// check interior rRESPA cutoff
if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3])
error->all(FLERR,"Pair cutoff < Respa interior cutoff");

2
src/CLASS2/pair_lj_class2_coul_long.h
View File

@ -40,7 +40,7 @@ class PairLJClass2CoulLong : public Pair {
void write_data(FILE *);
void write_data_all(FILE *);
double single(int, int, int, int, double, double, double, double &);
void compute_inner();
void compute_middle();
void compute_outer(int, int);

2
src/KOKKOS/comm_kokkos.cpp
View File

@ -210,7 +210,7 @@ void CommKokkos::forward_comm_device(int dummy)
MPI_Send(k_buf_send.view<DeviceType>().data(),
n,MPI_DOUBLE,sendproc[iswap],0,world);
}
if (size_forward_recv[iswap]) {
MPI_Wait(&request,MPI_STATUS_IGNORE);
atomKK->modified(ExecutionSpaceFromDevice<DeviceType>::

6
src/KOKKOS/domain_kokkos.cpp
View File

@ -340,11 +340,11 @@ struct DomainPBCFunctor {
void DomainKokkos::pbc()
{
if (lmp->kokkos->exchange_comm_classic) {
// reduce GPU data movement
atomKK->sync(Host,X_MASK|V_MASK|MASK_MASK|IMAGE_MASK);
Domain::pbc();
atomKK->modified(Host,X_MASK|V_MASK|MASK_MASK|IMAGE_MASK);

2
src/KOKKOS/kokkos.cpp
View File

@ -187,7 +187,7 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
binsize = 0.0;
#ifdef KOKKOS_ENABLE_CUDA
cuda_aware_flag = 1;
cuda_aware_flag = 1;
#else
cuda_aware_flag = 0;
#endif

2
src/KOKKOS/npair_kokkos.cpp
View File

@ -101,7 +101,7 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI,SIZE>::copy_stencil_info()
// copy stencil to device as it may have changed
int maxstencil = ns->get_maxstencil();
if (maxstencil > k_stencil.extent(0))
k_stencil = DAT::tdual_int_1d("neighlist:stencil",maxstencil);
for (int k = 0; k < maxstencil; k++)

4
src/KOKKOS/pair_kokkos.h
View File

@ -293,7 +293,7 @@ struct PairComputeFunctor {
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if(rsq < (STACKPARAMS?c.m_cutsq[itype][jtype]:c.d_cutsq(itype,jtype))) {
const F_FLOAT fpair = factor_lj*c.template compute_fpair<STACKPARAMS,Specialisation>(rsq,i,j,itype,jtype);
ftmp.x += delx*fpair;
@ -412,7 +412,7 @@ struct PairComputeFunctor {
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if(rsq < (STACKPARAMS?c.m_cutsq[itype][jtype]:c.d_cutsq(itype,jtype))) {
const F_FLOAT fpair = factor_lj*c.template compute_fpair<STACKPARAMS,Specialisation>(rsq,i,j,itype,jtype);
fev_tmp.f[0] += delx*fpair;

16
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -584,7 +584,7 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeDeidrj,const type
const int jj = team.league_rank() / ((inum+team.team_size()-1)/team.team_size());
const int ninside = d_ninside(ii);
if (jj >= ninside) return;
my_sna.compute_deidrj(team,ii,jj);
}
@ -619,9 +619,9 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeForce<NEIGHFLAG,E
a_f(j,0) -= fij[0];
a_f(j,1) -= fij[1];
a_f(j,2) -= fij[2];
// tally global and per-atom virial contribution
if (EVFLAG) {
if (vflag_either) {
v_tally_xyz<NEIGHFLAG>(ev,i,j,
@ -630,7 +630,7 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeForce<NEIGHFLAG,E
-my_sna.rij(ii,jj,2));
}
}
});
});
@ -649,16 +649,16 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeForce<NEIGHFLAG,E
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
double evdwl = d_coeffi[0];
// E = beta.B + 0.5*B^t.alpha.B
// linear contributions
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
evdwl += d_coeffi[icoeff+1]*my_sna.blist(ii,icoeff);
// quadratic contributions
if (quadraticflag) {
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {

2
src/KOKKOS/sna_kokkos.h
View File

@ -165,7 +165,7 @@ inline
void grow_rij(int, int);
int twojmax, diagonalstyle;
t_sna_2d blist;
t_sna_2c ulisttot;
t_sna_2c_lr ulisttot_lr;

30
src/KOKKOS/sna_kokkos_impl.h
View File

@ -83,7 +83,7 @@ void SNAKokkos<DeviceType>::build_indexlist()
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
h_idxcg_block(j1,j2,j) = idxcg_count;
h_idxcg_block(j1,j2,j) = idxcg_count;
for (int m1 = 0; m1 <= j1; m1++)
for (int m2 = 0; m2 <= j2; m2++)
idxcg_count++;
@ -98,9 +98,9 @@ void SNAKokkos<DeviceType>::build_indexlist()
auto h_idxu_block = Kokkos::create_mirror_view(idxu_block);
int idxu_count = 0;
for(int j = 0; j <= twojmax; j++) {
h_idxu_block[j] = idxu_count;
h_idxu_block[j] = idxu_count;
for(int mb = 0; mb <= j; mb++)
for(int ma = 0; ma <= j; ma++)
idxu_count++;
@ -110,16 +110,16 @@ void SNAKokkos<DeviceType>::build_indexlist()
// index list for beta and B
int idxb_count = 0;
int idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) idxb_count++;
idxb_max = idxb_count;
idxb = Kokkos::View<int*[3], DeviceType>("SNAKokkos::idxb",idxb_max);
auto h_idxb = Kokkos::create_mirror_view(idxb);
idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
@ -142,7 +142,7 @@ void SNAKokkos<DeviceType>::build_indexlist()
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
if (j >= j1) {
h_idxb_block(j1,j2,j) = idxb_count;
h_idxb_block(j1,j2,j) = idxb_count;
idxb_count++;
}
}
@ -158,19 +158,19 @@ void SNAKokkos<DeviceType>::build_indexlist()
for (int mb = 0; 2*mb <= j; mb++)
for (int ma = 0; ma <= j; ma++)
idxz_count++;
idxz_max = idxz_count;
idxz = Kokkos::View<int*[10], DeviceType>("SNAKokkos::idxz",idxz_max);
auto h_idxz = Kokkos::create_mirror_view(idxz);
idxz_block = Kokkos::View<int***, DeviceType>("SNAKokkos::idxz_block", jdim,jdim,jdim);
auto h_idxz_block = Kokkos::create_mirror_view(idxz_block);
idxz_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
h_idxz_block(j1,j2,j) = idxz_count;
h_idxz_block(j1,j2,j) = idxz_count;
// find right beta(ii,jjb) entry
// multiply and divide by j+1 factors
@ -226,7 +226,7 @@ void SNAKokkos<DeviceType>::grow_rij(int newnatom, int newnmax)
blist = t_sna_2d("sna:blist",natom,idxb_max);
ulisttot = t_sna_2c("sna:ulisttot",natom,idxu_max);
if (!Kokkos::Impl::is_same<typename DeviceType::array_layout,Kokkos::LayoutRight>::value)
if (!Kokkos::Impl::is_same<typename DeviceType::array_layout,Kokkos::LayoutRight>::value)
ulisttot_lr = t_sna_2c_lr("sna:ulisttot_lr",natom,idxu_max);
zlist = t_sna_2c("sna:zlist",natom,idxz_max);
@ -306,7 +306,7 @@ void SNAKokkos<DeviceType>::compute_zi(const int& iter)
const double* cgblock = cglist.data() + idxcg_block(j1,j2,j);
zlist(iatom,jjz).re = 0.0;
zlist(iatom,jjz).re = 0.0;
zlist(iatom,jjz).im = 0.0;
int jju1 = idxu_block[j1] + (j1+1)*mb1min;
@ -419,7 +419,7 @@ void SNAKokkos<DeviceType>::compute_yi(int iter,
if (j1 == j) {
if (j2 == j) betaj = 3*beta(iatom,jjb);
else betaj = 2*beta(iatom,jjb);
} else betaj = beta(iatom,jjb);
} else betaj = beta(iatom,jjb);
} else if (j >= j2) {
const int jjb = idxb_block(j,j2,j1);
if (j2 == j) betaj = 2*beta(iatom,jjb)*(j1+1)/(j+1.0);
@ -1176,7 +1176,7 @@ void SNAKokkos<DeviceType>::init_clebsch_gordan()
factorial((j + cc2) / 2) *
factorial((j - cc2) / 2) *
(j + 1));
h_cglist[idxcg_count] = sum * dcg * sfaccg;
idxcg_count++;
}
@ -1278,7 +1278,7 @@ double SNAKokkos<DeviceType>::memory_usage()
if (!Kokkos::Impl::is_same<typename DeviceType::array_layout,Kokkos::LayoutRight>::value)
bytes += natom * idxu_max * sizeof(double) * 2; // ulisttot_lr
bytes += natom * idxu_max * 3 * sizeof(double) * 2; // dulist
bytes += natom * idxz_max * sizeof(double) * 2; // zlist
bytes += natom * idxb_max * sizeof(double); // blist
bytes += natom * idxu_max * sizeof(double) * 2; // ylist

22
src/KSPACE/ewald_dipole.cpp
View File

@ -167,7 +167,7 @@ void EwaldDipole::init()
NewtonSolve(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
if (g_ewald_new > 0.0) g_ewald = g_ewald_new;
else error->warning(FLERR,"Ewald/disp Newton solver failed, "
"using old method to estimate g_ewald");
"using old method to estimate g_ewald");
}
// setup EwaldDipole coefficients so can print stats
@ -246,7 +246,7 @@ void EwaldDipole::setup()
double err;
kxmax = 1;
kymax = 1;
kzmax = 1;
kzmax = 1;
// set kmax in 3 directions to respect accuracy
@ -462,7 +462,7 @@ void EwaldDipole::compute(int eflag, int vflag)
vc[k][4] += vcik[4] = -(partial_peratom * mu[i][0] * eg[k][2]);
vc[k][5] += vcik[5] = -(partial_peratom * mu[i][1] * eg[k][2]);
// taking re-part of struct_fact x exp(i*k*ri)
// taking re-part of struct_fact x exp(i*k*ri)
// (for per-atom energy and virial calc.)
if (evflag_atom) {
@ -653,12 +653,12 @@ void EwaldDipole::eik_dot_r()
muz = mu[i][2];
// dir 1: (0,l,m)
mudotk = (muy*l*unitk[1] + muz*m*unitk[2]);
mudotk = (muy*l*unitk[1] + muz*m*unitk[2]);
cstr1 += mudotk*(cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i]);
sstr1 += mudotk*(sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i]);
// dir 2: (0,l,-m)
mudotk = (muy*l*unitk[1] - muz*m*unitk[2]);
mudotk = (muy*l*unitk[1] - muz*m*unitk[2]);
cstr2 += mudotk*(cs[l][1][i]*cs[m][2][i]+sn[l][1][i]*sn[m][2][i]);
sstr2 += mudotk*(sn[l][1][i]*cs[m][2][i]-cs[l][1][i]*sn[m][2][i]);
}
@ -685,12 +685,12 @@ void EwaldDipole::eik_dot_r()
muz = mu[i][2];
// dir 1: (k,0,m)
mudotk = (mux*k*unitk[0] + muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] + muz*m*unitk[2]);
cstr1 += mudotk*(cs[k][0][i]*cs[m][2][i]-sn[k][0][i]*sn[m][2][i]);
sstr1 += mudotk*(sn[k][0][i]*cs[m][2][i]+cs[k][0][i]*sn[m][2][i]);
// dir 2: (k,0,-m)
mudotk = (mux*k*unitk[0] - muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] - muz*m*unitk[2]);
cstr2 += mudotk*(cs[k][0][i]*cs[m][2][i]+sn[k][0][i]*sn[m][2][i]);
sstr2 += mudotk*(sn[k][0][i]*cs[m][2][i]-cs[k][0][i]*sn[m][2][i]);
}
@ -724,28 +724,28 @@ void EwaldDipole::eik_dot_r()
muz = mu[i][2];
// dir 1: (k,l,m)
mudotk = (mux*k*unitk[0] + muy*l*unitk[1] + muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] + muy*l*unitk[1] + muz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i];
slpm = sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i];
cstr1 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr1 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 2: (k,-l,m)
mudotk = (mux*k*unitk[0] - muy*l*unitk[1] + muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] - muy*l*unitk[1] + muz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i];
slpm = -sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i];
cstr2 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr2 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 3: (k,l,-m)
mudotk = (mux*k*unitk[0] + muy*l*unitk[1] - muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] + muy*l*unitk[1] - muz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i];
slpm = sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i];
cstr3 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr3 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 4: (k,-l,-m)
mudotk = (mux*k*unitk[0] - muy*l*unitk[1] - muz*m*unitk[2]);
mudotk = (mux*k*unitk[0] - muy*l*unitk[1] - muz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i];
slpm = -sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i];
cstr4 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);

26
src/KSPACE/ewald_dipole_spin.cpp
View File

@ -36,7 +36,7 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
EwaldDipoleSpin::EwaldDipoleSpin(LAMMPS *lmp) :
EwaldDipoleSpin::EwaldDipoleSpin(LAMMPS *lmp) :
EwaldDipole(lmp)
{
dipoleflag = 0;
@ -157,7 +157,7 @@ void EwaldDipoleSpin::init()
NewtonSolve(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
if (g_ewald_new > 0.0) g_ewald = g_ewald_new;
else error->warning(FLERR,"Ewald/disp Newton solver failed, "
"using old method to estimate g_ewald");
"using old method to estimate g_ewald");
}
// setup EwaldDipoleSpin coefficients so can print stats
@ -236,7 +236,7 @@ void EwaldDipoleSpin::setup()
double err;
kxmax = 1;
kymax = 1;
kzmax = 1;
kzmax = 1;
// set kmax in 3 directions to respect accuracy
@ -440,7 +440,7 @@ void EwaldDipoleSpin::compute(int eflag, int vflag)
vc[k][4] += vcik[4] = -(partial_peratom * spx * eg[k][2]);
vc[k][5] += vcik[5] = -(partial_peratom * spy * eg[k][2]);
// taking re-part of struct_fact x exp(i*k*ri)
// taking re-part of struct_fact x exp(i*k*ri)
// (for per-atom energy and virial calc.)
if (evflag_atom) {
@ -639,12 +639,12 @@ void EwaldDipoleSpin::eik_dot_r()
spz = sp[i][2]*sp[i][3];
// dir 1: (0,l,m)
mudotk = (spy*l*unitk[1] + spz*m*unitk[2]);
mudotk = (spy*l*unitk[1] + spz*m*unitk[2]);
cstr1 += mudotk*(cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i]);
sstr1 += mudotk*(sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i]);
// dir 2: (0,l,-m)
mudotk = (spy*l*unitk[1] - spz*m*unitk[2]);
mudotk = (spy*l*unitk[1] - spz*m*unitk[2]);
cstr2 += mudotk*(cs[l][1][i]*cs[m][2][i]+sn[l][1][i]*sn[m][2][i]);
sstr2 += mudotk*(sn[l][1][i]*cs[m][2][i]-cs[l][1][i]*sn[m][2][i]);
}
@ -671,12 +671,12 @@ void EwaldDipoleSpin::eik_dot_r()
spz = sp[i][2]*sp[i][3];
// dir 1: (k,0,m)
mudotk = (spx*k*unitk[0] + spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] + spz*m*unitk[2]);
cstr1 += mudotk*(cs[k][0][i]*cs[m][2][i]-sn[k][0][i]*sn[m][2][i]);
sstr1 += mudotk*(sn[k][0][i]*cs[m][2][i]+cs[k][0][i]*sn[m][2][i]);
// dir 2: (k,0,-m)
mudotk = (spx*k*unitk[0] - spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] - spz*m*unitk[2]);
cstr2 += mudotk*(cs[k][0][i]*cs[m][2][i]+sn[k][0][i]*sn[m][2][i]);
sstr2 += mudotk*(sn[k][0][i]*cs[m][2][i]-cs[k][0][i]*sn[m][2][i]);
}
@ -710,28 +710,28 @@ void EwaldDipoleSpin::eik_dot_r()
spz = sp[i][2]*sp[i][3];
// dir 1: (k,l,m)
mudotk = (spx*k*unitk[0] + spy*l*unitk[1] + spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] + spy*l*unitk[1] + spz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i];
slpm = sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i];
cstr1 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr1 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 2: (k,-l,m)
mudotk = (spx*k*unitk[0] - spy*l*unitk[1] + spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] - spy*l*unitk[1] + spz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i];
slpm = -sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i];
cstr2 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr2 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 3: (k,l,-m)
mudotk = (spx*k*unitk[0] + spy*l*unitk[1] - spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] + spy*l*unitk[1] - spz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i];
slpm = sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i];
cstr3 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
sstr3 += mudotk*(sn[k][0][i]*clpm + cs[k][0][i]*slpm);
// dir 4: (k,-l,-m)
mudotk = (spx*k*unitk[0] - spy*l*unitk[1] - spz*m*unitk[2]);
mudotk = (spx*k*unitk[0] - spy*l*unitk[1] - spz*m*unitk[2]);
clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i];
slpm = -sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i];
cstr4 += mudotk*(cs[k][0][i]*clpm - sn[k][0][i]*slpm);
@ -768,7 +768,7 @@ void EwaldDipoleSpin::slabcorr()
double spz;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
spz = sp[i][2]*sp[i][3];
spin += spz;
}

6
src/KSPACE/ewald_dipole_spin.h
View File

@ -33,13 +33,13 @@ class EwaldDipoleSpin : public EwaldDipole {
void compute(int, int);
protected:
double hbar; // reduced Planck's constant
double mub; // Bohr's magneton
double hbar; // reduced Planck's constant
double mub; // Bohr's magneton
double mu_0; // vacuum permeability
double mub2mu0; // prefactor for mech force
double mub2mu0hbinv; // prefactor for mag force
void spsum_musq();
void spsum_musq();
virtual void eik_dot_r();
void slabcorr();

16
src/KSPACE/pppm_dipole.cpp
View File

@ -58,19 +58,19 @@ enum{FORWARD_MU,FORWARD_MU_PERATOM};
/* ---------------------------------------------------------------------- */
PPPMDipole::PPPMDipole(LAMMPS *lmp) : PPPM(lmp),
densityx_brick_dipole(NULL), densityy_brick_dipole(NULL),
densityx_brick_dipole(NULL), densityy_brick_dipole(NULL),
densityz_brick_dipole(NULL),
vdxx_brick_dipole(NULL), vdyy_brick_dipole(NULL), vdzz_brick_dipole(NULL),
vdxy_brick_dipole(NULL), vdxz_brick_dipole(NULL), vdyz_brick_dipole(NULL),
ux_brick_dipole(NULL), uy_brick_dipole(NULL), uz_brick_dipole(NULL),
v0x_brick_dipole(NULL), v1x_brick_dipole(NULL),
v2x_brick_dipole(NULL), v3x_brick_dipole(NULL), v4x_brick_dipole(NULL),
v5x_brick_dipole(NULL), v0y_brick_dipole(NULL), v1y_brick_dipole(NULL),
v2y_brick_dipole(NULL), v3y_brick_dipole(NULL), v4y_brick_dipole(NULL),
v5y_brick_dipole(NULL), v0z_brick_dipole(NULL), v1z_brick_dipole(NULL),
v2z_brick_dipole(NULL), v3z_brick_dipole(NULL), v4z_brick_dipole(NULL),
v5z_brick_dipole(NULL), work3(NULL), work4(NULL),
densityx_fft_dipole(NULL), densityy_fft_dipole(NULL),
v2x_brick_dipole(NULL), v3x_brick_dipole(NULL), v4x_brick_dipole(NULL),
v5x_brick_dipole(NULL), v0y_brick_dipole(NULL), v1y_brick_dipole(NULL),
v2y_brick_dipole(NULL), v3y_brick_dipole(NULL), v4y_brick_dipole(NULL),
v5y_brick_dipole(NULL), v0z_brick_dipole(NULL), v1z_brick_dipole(NULL),
v2z_brick_dipole(NULL), v3z_brick_dipole(NULL), v4z_brick_dipole(NULL),
v5z_brick_dipole(NULL), work3(NULL), work4(NULL),
densityx_fft_dipole(NULL), densityy_fft_dipole(NULL),
densityz_fft_dipole(NULL)
{
dipoleflag = 1;

2
src/KSPACE/pppm_dipole.h
View File

@ -38,7 +38,7 @@ class PPPMDipole : public PPPM {
protected:
void set_grid_global();
double newton_raphson_f();
double newton_raphson_f();
void allocate();
void allocate_peratom();

8
src/KSPACE/pppm_dipole_spin.cpp
View File

@ -52,7 +52,7 @@ enum{FORWARD_MU,FORWARD_MU_PERATOM};
/* ---------------------------------------------------------------------- */
PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp) :
PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp) :
PPPMDipole(lmp)
{
dipoleflag = 0;
@ -147,7 +147,7 @@ void PPPMDipoleSpin::init()
// kspace TIP4P not yet supported
// qdist = offset only for TIP4P fictitious charge
qdist = 0.0;
qdist = 0.0;
if (tip4pflag)
error->all(FLERR,"Cannot yet use TIP4P with PPPMDipoleSpin");
@ -668,7 +668,7 @@ void PPPMDipoleSpin::slabcorr()
double spz;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
spz = sp[i][2]*sp[i][3];
spin += spz;
}
@ -729,7 +729,7 @@ void PPPMDipoleSpin::spsum_spsq()
spsqsum_local += spx*spx + spy*spy + spz*spz;
}
// store results into pppm_dipole quantities
// store results into pppm_dipole quantities
MPI_Allreduce(&spsum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&spsqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world);

4
src/KSPACE/pppm_dipole_spin.h
View File

@ -32,8 +32,8 @@ class PPPMDipoleSpin : public PPPMDipole {
void compute(int, int);
protected:
double hbar; // reduced Planck's constant
double mub; // Bohr's magneton
double hbar; // reduced Planck's constant
double mub; // Bohr's magneton
double mu_0; // vacuum permeability
double mub2mu0; // prefactor for mech force
double mub2mu0hbinv; // prefactor for mag force

8
src/MANYBODY/pair_airebo.cpp
View File

@ -3638,9 +3638,9 @@ void PairAIREBO::read_file(char *filename)
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
if (1 != sscanf(s,"%lg",&reqM_HH)) ++cerror;
}
}
// check for errors parsing global parameters
MPI_Bcast(&cerror,1,MPI_INT,0,world);
@ -3654,7 +3654,7 @@ void PairAIREBO::read_file(char *filename)
cerror = numpar = 0;
if (me == 0) {
// gC spline
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
@ -3899,7 +3899,7 @@ void PairAIREBO::read_file(char *filename)
fclose(fp);
}
// check for errors parsing spline data
MPI_Bcast(&cerror,1,MPI_INT,0,world);

2
src/RIGID/rigid_const.h
View File

@ -32,7 +32,7 @@ namespace LAMMPS_NS {
ANGMOM = 1<<7,
TORQUE = 1<<8
};
static const double TOLERANCE = 1.0e-6;
static const double EPSILON = 1.0e-7;
static const double BIG = 1.0e20;

4
src/SNAP/pair_snap.cpp
View File

@ -108,7 +108,7 @@ void PairSNAP::compute(int eflag, int vflag)
// compute dE_i/dB_i = beta_i for all i in list
if (quadraticflag || eflag)
if (quadraticflag || eflag)
compute_bispectrum();
compute_beta();
@ -165,7 +165,7 @@ void PairSNAP::compute(int eflag, int vflag)
snaptr->compute_ui(ninside);
// for neighbors of I within cutoff:
// compute Fij = dEi/dRj = -dEi/dRi
// compute Fij = dEi/dRj = -dEi/dRi
// add to Fi, subtract from Fj
snaptr->compute_yi(beta[ii]);

40
src/SNAP/sna.cpp
View File

@ -171,7 +171,7 @@ void SNA::build_indexlist()
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
idxcg_block[j1][j2][j] = idxcg_count;
idxcg_block[j1][j2][j] = idxcg_count;
for (int m1 = 0; m1 <= j1; m1++)
for (int m2 = 0; m2 <= j2; m2++)
idxcg_count++;
@ -185,9 +185,9 @@ void SNA::build_indexlist()
"sna:idxu_block");
int idxu_count = 0;
for(int j = 0; j <= twojmax; j++) {
idxu_block[j] = idxu_count;
idxu_block[j] = idxu_count;
for(int mb = 0; mb <= j; mb++)
for(int ma = 0; ma <= j; ma++)
idxu_count++;
@ -196,15 +196,15 @@ void SNA::build_indexlist()
// index list for beta and B
int idxb_count = 0;
int idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) idxb_count++;
idxb_max = idxb_count;
idxb = new SNA_BINDICES[idxb_max];
idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
@ -225,7 +225,7 @@ void SNA::build_indexlist()
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
if (j >= j1) {
idxb_block[j1][j2][j] = idxb_count;
idxb_block[j1][j2][j] = idxb_count;
idxb_count++;
}
}
@ -240,18 +240,18 @@ void SNA::build_indexlist()
for (int mb = 0; 2*mb <= j; mb++)
for (int ma = 0; ma <= j; ma++)
idxz_count++;
idxz_max = idxz_count;
idxz = new SNA_ZINDICES[idxz_max];
memory->create(idxz_block, jdim, jdim, jdim,
"sna:idxz_block");
idxz_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
idxz_block[j1][j2][j] = idxz_count;
idxz_block[j1][j2][j] = idxz_count;
// find right beta[jjb] entry
// multiply and divide by j+1 factors
@ -481,7 +481,7 @@ void SNA::compute_yi(const double* beta)
if (j1 == j) {
if (j2 == j) betaj = 3*beta[jjb];
else betaj = 2*beta[jjb];
} else betaj = beta[jjb];
} else betaj = beta[jjb];
} else if (j >= j2) {
const int jjb = idxb_block[j][j2][j1];
if (j2 == j) betaj = 2*beta[jjb]*(j1+1)/(j+1.0);
@ -549,7 +549,7 @@ void SNA::compute_deidrj(double* dedr)
double jjjmambyarray_i = ylist_i[jju];
for(int k = 0; k < 3; k++)
dedr[k] +=
dedr[k] +=
(dudr_r[k] * jjjmambyarray_r +
dudr_i[k] * jjjmambyarray_i)*0.5;
jju++;
@ -588,24 +588,24 @@ void SNA::compute_bi()
double sumzu = 0.0;
for (int mb = 0; 2*mb < j; mb++)
for (int ma = 0; ma <= j; ma++) {
sumzu += ulisttot_r[jju]*zlist_r[jjz] +
sumzu += ulisttot_r[jju]*zlist_r[jjz] +
ulisttot_i[jju]*zlist_i[jjz];
jjz++;
jju++;
} // end loop over ma, mb
} // end loop over ma, mb
// For j even, handle middle column
if (j%2 == 0) {
int mb = j/2;
for(int ma = 0; ma < mb; ma++) {
sumzu += ulisttot_r[jju]*zlist_r[jjz] +
sumzu += ulisttot_r[jju]*zlist_r[jjz] +
ulisttot_i[jju]*zlist_i[jjz];
jjz++;
jju++;
}
sumzu += 0.5*(ulisttot_r[jju]*zlist_r[jjz] +
sumzu += 0.5*(ulisttot_r[jju]*zlist_r[jjz] +
ulisttot_i[jju]*zlist_i[jjz]);
} // end if jeven
@ -1485,7 +1485,7 @@ void SNA::init_clebsch_gordan()
factorial((j - j2 + aa2) / 2 + z) *
factorial((j - j1 - bb2) / 2 + z));
}
cc2 = 2 * m - j;
dcg = deltacg(j1, j2, j);
sfaccg = sqrt(factorial((j1 + aa2) / 2) *
@ -1495,7 +1495,7 @@ void SNA::init_clebsch_gordan()
factorial((j + cc2) / 2) *
factorial((j - cc2) / 2) *
(j + 1));
cglist[idxcg_count] = sum * dcg * sfaccg;
idxcg_count++;
}
@ -1519,7 +1519,7 @@ void SNA::print_clebsch_gordan()
for (int j1 = 0; j1 <= twojmax; j1++)
for (int j2 = 0; j2 <= j1; j2++)
if (j1-j2 <= j && j1+j2 >= j && (j1+j2+j)%2 == 0) {
int idxcg_count = idxcg_block[j1][j2][j];
int idxcg_count = idxcg_block[j1][j2][j];
for (int m1 = 0; m1 <= j1; m1++) {
aa2 = 2*m1-j1;
for (int m2 = 0; m2 <= j2; m2++) {

4
src/SNAP/sna.h
View File

@ -81,8 +81,8 @@ private:
int idxcg_max, idxu_max, idxz_max, idxb_max;
double** rootpqarray;
double* cglist;
int*** idxcg_block;
double* cglist;
int*** idxcg_block;
double* ulisttot_r, * ulisttot_i;
double** ulist_r_ij, ** ulist_i_ij;

2
src/SPIN/atom_vec_spin.h
View File

@ -68,7 +68,7 @@ class AtomVecSpin : public AtomVec {
int *type,*mask;
imageint *image;
double **x,**v,**f; // lattice quantities
// spin quantities
double **sp; // sp[i][0-2] direction of the spin i
// sp[i][3] atomic magnetic moment of the spin i

2
src/SPIN/fix_nve_spin.cpp
View File

@ -247,7 +247,7 @@ void FixNVESpin::init()
locksetforcespin = (FixSetForceSpin *) modify->fix[iforce];
}
}
// setting the sector variables/lists
nsectors = 0;

10
src/SPIN/fix_precession_spin.cpp
View File

@ -126,7 +126,7 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
nay *= inorm;
naz *= inorm;
}
if (cubic_flag) {
inorm = 1.0/sqrt(nc1x*nc1x + nc1y*nc1y + nc1z*nc1z);
nc1x *= inorm;
@ -244,7 +244,7 @@ void FixPrecessionSpin::post_force(int /* vflag */)
double **sp = atom->sp;
const int nlocal = atom->nlocal;
double spi[3], fmi[3], epreci;
eflag = 0;
eprec = 0.0;
for (int i = 0; i < nlocal; i++) {
@ -317,7 +317,7 @@ double FixPrecessionSpin::compute_anisotropy_energy(double spi[3])
double energy = 0.0;
double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
energy = Ka*scalar*scalar;
return energy;
return energy;
}
/* ---------------------------------------------------------------------- */
@ -391,11 +391,11 @@ void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3])
six1 = 2.0*skx*sky2*skz2;
six2 = 2.0*sky*skx2*skz2;
six3 = 2.0*skz*skx2*sky2;
sixx = k2ch*(nc1x*six1 + nc2x*six2 + nc3x*six3);
sixy = k2ch*(nc1y*six1 + nc2y*six2 + nc3y*six3);
sixz = k2ch*(nc1z*six1 + nc2z*six2 + nc3z*six3);
fmi[0] += fourx + sixx;
fmi[1] += foury + sixy;
fmi[2] += fourz + sixz;

4
src/SPIN/fix_precession_spin.h
View File

@ -42,7 +42,7 @@ class FixPrecessionSpin : public Fix {
int zeeman_flag, aniso_flag, cubic_flag;
void compute_single_precession(int, double *, double *);
void compute_zeeman(int, double *);
// uniaxial aniso calculations
void compute_anisotropy(double *, double *);
@ -52,7 +52,7 @@ class FixPrecessionSpin : public Fix {
void compute_cubic(double *, double *);
double compute_cubic_energy(double *);
protected:
int style; // style of the magnetic precession

2
src/SPIN/fix_setforce_spin.cpp
View File

@ -140,7 +140,7 @@ void FixSetForceSpin::single_setforce_spin(int i, double fmi[3])
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
force_flag = 0;
// constant force
if (varflag == CONSTANT) {

2
src/SPIN/fix_setforce_spin.h
View File

@ -29,7 +29,7 @@ class FixSetForceSpin : public FixSetForce {
FixSetForceSpin(class LAMMPS *, int, char **);
virtual void post_force(int);
void post_force_respa(int, int, int);
void single_setforce_spin(int, double *);
void single_setforce_spin(int, double *);
};
}

26
src/SPIN/min_spin_cg.cpp
View File

@ -16,8 +16,8 @@
Julien Tranchida (SNL)
Please cite the related publication:
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
preprint arXiv:1904.02669.
------------------------------------------------------------------------- */
@ -105,7 +105,7 @@ void MinSpinCG::init()
error->warning(FLERR,"Line search incompatible gneb");
// set back use_line_search to 0 if more than one replica
if (linestyle == 3 && nreplica == 1){
use_line_search = 1;
}
@ -201,10 +201,10 @@ int MinSpinCG::iterate(int maxiter)
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
@ -249,7 +249,7 @@ int MinSpinCG::iterate(int maxiter)
// energy tolerance criterion
// only check after DELAYSTEP elapsed since velocties reset to 0
// sync across replicas if running multi-replica minimization
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
@ -365,7 +365,7 @@ void MinSpinCG::calc_search_direction()
MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,world);
// Sum over all replicas. Good for GNEB.
if (nreplica > 1) {
g2 = g2_global * factor;
g2old = g2old_global * factor;
@ -374,9 +374,9 @@ void MinSpinCG::calc_search_direction()
}
if (fabs(g2_global) < 1.0e-60) beta = 0.0;
else beta = g2_global / g2old_global;
// calculate conjugate direction
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = (beta * p_s[i] - g_cur[i]) * factor;
g_old[i] = g_cur[i] * factor;
@ -401,9 +401,9 @@ void MinSpinCG::advance_spins()
for (int i = 0; i < nlocal; i++) {
rodrigues_rotation(p_s + 3 * i, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
@ -414,7 +414,7 @@ void MinSpinCG::advance_spins()
(R. Murray, Z. Li, and S. Shankar Sastry,
A Mathematical Introduction to
Robotic Manipulation (1994), p. 28 and 30).
upp_tr - vector x, y, z so that one calculate
U = exp(A) with A= [[0, x, y],
[-x, 0, z],
@ -431,7 +431,7 @@ void MinSpinCG::rodrigues_rotation(const double *upp_tr, double *out)
fabs(upp_tr[2]) < 1.0e-40){
// if upp_tr is zero, return unity matrix
for(int k = 0; k < 3; k++){
for(int m = 0; m < 3; m++){
if (m == k) out[3 * k + m] = 1.0;

16
src/SPIN/min_spin_lbfgs.cpp
View File

@ -16,8 +16,8 @@
Julien Tranchida (SNL)
Please cite the related publication:
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
preprint arXiv:1904.02669.
------------------------------------------------------------------------- */
@ -213,10 +213,10 @@ int MinSpinLBFGS::iterate(int maxiter)
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
@ -264,7 +264,7 @@ int MinSpinLBFGS::iterate(int maxiter)
// energy tolerance criterion
// only check after DELAYSTEP elapsed since velocties reset to 0
// sync across replicas if running multi-replica minimization
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
@ -526,9 +526,9 @@ void MinSpinLBFGS::advance_spins()
for (int i = 0; i < nlocal; i++) {
rodrigues_rotation(p_s + 3 * i, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
@ -539,7 +539,7 @@ void MinSpinLBFGS::advance_spins()
(R. Murray, Z. Li, and S. Shankar Sastry,
A Mathematical Introduction to
Robotic Manipulation (1994), p. 28 and 30).
upp_tr - vector x, y, z so that one calculate
U = exp(A) with A= [[0, x, y],
[-x, 0, z],

2
src/SPIN/pair_spin.cpp
View File

@ -82,7 +82,7 @@ void PairSpin::init_style()
bool have_fix = ((modify->find_fix_by_style("^nve/spin") != -1)
|| (modify->find_fix_by_style("^neb/spin") != -1));
if (!have_fix && (comm->me == 0))
error->warning(FLERR,"Using spin pair style without nve/spin or neb/spin");

88
src/SPIN/pair_spin_dipole_cut.cpp
View File

@ -76,9 +76,9 @@ void PairSpinDipoleCut::settings(int narg, char **arg)
PairSpin::settings(narg,arg);
cut_spin_long_global = force->numeric(FLERR,arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
@ -99,10 +99,10 @@ void PairSpinDipoleCut::settings(int narg, char **arg)
void PairSpinDipoleCut::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (narg != 3)
if (narg != 3)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
@ -128,9 +128,9 @@ void PairSpinDipoleCut::coeff(int narg, char **arg)
double PairSpinDipoleCut::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cut_spin_long[j][i] = cut_spin_long[i][j];
return cut_spin_long_global;
}
@ -163,8 +163,8 @@ void *PairSpinDipoleCut::extract(const char *str, int &dim)
void PairSpinDipoleCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rinv,r2inv,r3inv,rsq,local_cut2,evdwl,ecoul;
double xi[3],rij[3],eij[3],spi[4],spj[4],fi[3],fmi[3];
@ -172,13 +172,13 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
int *type = atom->type;
int nlocal = atom->nlocal;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
double **sp = atom->sp;
inum = list->inum;
ilist = list->ilist;
@ -194,9 +194,9 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
jnum = numneigh[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
spi[3] = sp[i][3];
itype = type[i];
@ -206,15 +206,15 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
evdwl = 0.0;
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
@ -229,23 +229,23 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r3inv = r2inv*rinv;
compute_dipolar(i,j,eij,fmi,spi,spj,r3inv);
if (lattice_flag) compute_dipolar_mech(i,j,eij,fi,spi,spj,r2inv);
}
// force accumulation
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
@ -269,21 +269,21 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
{
int j,jnum,itype,jtype,ntypes;
int *jlist,*numneigh,**firstneigh;
int j,jnum,itype,jtype,ntypes;
int *jlist,*numneigh,**firstneigh;
double rsq,rinv,r2inv,r3inv,local_cut2;
double xi[3],rij[3],eij[3],spi[4],spj[4];
int k,locflag;
int *type = atom->type;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double **sp = atom->sp;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
ntypes = atom->ntypes;
locflag = 0;
@ -307,28 +307,28 @@ void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[2] = sp[ii][2];
spi[3] = sp[ii][3];
jlist = firstneigh[ii];
jnum = numneigh[ii];
jnum = numneigh[ii];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
@ -344,9 +344,9 @@ void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r3inv = r2inv*rinv;
// compute dipolar interaction
compute_dipolar(ii,j,eij,fmi,spi,spj,r3inv);
}
}
@ -357,7 +357,7 @@ void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
compute dipolar interaction between spins i and j
------------------------------------------------------------------------- */
void PairSpinDipoleCut::compute_dipolar(int /* i */, int /* j */, double eij[3],
void PairSpinDipoleCut::compute_dipolar(int /* i */, int /* j */, double eij[3],
double fmi[3], double spi[4], double spj[4], double r3inv)
{
double sjdotr;
@ -373,7 +373,7 @@ void PairSpinDipoleCut::compute_dipolar(int /* i */, int /* j */, double eij[3],
}
/* ----------------------------------------------------------------------
compute the mechanical force due to the dipolar interaction between
compute the mechanical force due to the dipolar interaction between
atom i and atom j
------------------------------------------------------------------------- */
@ -387,7 +387,7 @@ void PairSpinDipoleCut::compute_dipolar_mech(int /* i */, int /* j */, double ei
sisj = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2];
sieij = spi[0]*eij[0] + spi[1]*eij[1] + spi[2]*eij[2];
sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2];
bij = sisj - 5.0*sieij*sjeij;
pre = 3.0*mub2mu0*gigjri4;

12
src/SPIN/pair_spin_dipole_cut.h
View File

@ -34,22 +34,22 @@ class PairSpinDipoleCut : public PairSpin {
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void *extract(const char *, int &);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_dipolar(int, int, double *, double *, double *,
void compute_dipolar(int, int, double *, double *, double *,
double *, double);
void compute_dipolar_mech(int, int, double *, double *, double *,
void compute_dipolar_mech(int, int, double *, double *, double *,
double *, double);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double cut_spin_long_global; // global long cutoff distance
double cut_spin_long_global; // global long cutoff distance
protected:
double hbar; // reduced Planck's constant

86
src/SPIN/pair_spin_dipole_long.cpp
View File

@ -81,9 +81,9 @@ void PairSpinDipoleLong::settings(int narg, char **arg)
PairSpin::settings(narg,arg);
cut_spin_long_global = force->numeric(FLERR,arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
@ -103,10 +103,10 @@ void PairSpinDipoleLong::settings(int narg, char **arg)
void PairSpinDipoleLong::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (narg != 3)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
@ -148,9 +148,9 @@ void PairSpinDipoleLong::init_style()
double PairSpinDipoleLong::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cut_spin_long[j][i] = cut_spin_long[i][j];
return cut_spin_long_global;
}
@ -183,7 +183,7 @@ void *PairSpinDipoleLong::extract(const char *str, int &dim)
void PairSpinDipoleLong::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int i,j,ii,jj,inum,jnum,itype,jtype;
double r,rinv,r2inv,rsq;
double grij,expm2,t,erfc;
double evdwl,ecoul;
@ -193,7 +193,7 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
double fi[3],fmi[3];
double local_cut2;
double pre1,pre2,pre3;
int *ilist,*jlist,*numneigh,**firstneigh;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -202,9 +202,9 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
@ -225,9 +225,9 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
jnum = numneigh[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
spi[3] = sp[i][3];
itype = type[i];
@ -237,17 +237,17 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
evdwl = 0.0;
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
bij[0] = bij[1] = bij[2] = bij[3] = 0.0;
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
@ -279,22 +279,22 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
// force accumulation
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
if (eflag) {
if (rsq <= local_cut2) {
evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] +
evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] +
spi[2]*fmi[2];
evdwl *= hbar;
}
@ -314,21 +314,21 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
{
int j,jj,jnum,itype,jtype,ntypes;
int j,jj,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
int *jlist,*numneigh,**firstneigh;
double r,rinv,r2inv,rsq,grij,expm2,t,erfc;
double local_cut2,pre1,pre2,pre3;
double bij[4],xi[3],rij[3],eij[3],spi[4],spj[4];
int *type = atom->type;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double **sp = atom->sp;
double **fm_long = atom->fm_long;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
@ -362,16 +362,16 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
// computation of the exchange interaction
// loop over neighbors of atom i
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[2] = sp[ii][2];
spi[3] = sp[ii][3];
jlist = firstneigh[ii];
jnum = numneigh[ii];
jnum = numneigh[ii];
//itype = type[i];
for (jj = 0; jj < jnum; jj++) {
@ -379,14 +379,14 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
fmi[0] = fmi[1] = fmi[2] = 0.0;
bij[0] = bij[1] = bij[2] = bij[3] = 0.0;
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
@ -417,7 +417,7 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
}
// adding the kspace components to fm
fmi[0] += fm_long[ii][0];
fmi[1] += fm_long[ii][1];
fmi[2] += fm_long[ii][2];
@ -428,7 +428,7 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
compute dipolar interaction between spins i and j
------------------------------------------------------------------------- */
void PairSpinDipoleLong::compute_long(int /* i */, int /* j */, double eij[3],
void PairSpinDipoleLong::compute_long(int /* i */, int /* j */, double eij[3],
double bij[4], double fmi[3], double spi[4], double spj[4])
{
double sjeij,pre;
@ -447,7 +447,7 @@ void PairSpinDipoleLong::compute_long(int /* i */, int /* j */, double eij[3],
}
/* ----------------------------------------------------------------------
compute the mechanical force due to the dipolar interaction between
compute the mechanical force due to the dipolar interaction between
atom i and atom j
------------------------------------------------------------------------- */

12
src/SPIN/pair_spin_dipole_long.h
View File

@ -35,22 +35,22 @@ class PairSpinDipoleLong : public PairSpin {
void coeff(int, char **);
void init_style();
double init_one(int, int);
void *extract(const char *, int &);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_long(int, int, double *, double *, double *,
void compute_long(int, int, double *, double *, double *,
double *, double *);
void compute_long_mech(int, int, double *, double *, double *,
void compute_long_mech(int, int, double *, double *, double *,
double *, double *);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double cut_spin_long_global; // global long cutoff distance
double cut_spin_long_global; // global long cutoff distance
protected:
double hbar; // reduced Planck's constant

2
src/USER-CGDNA/pair_oxdna2_excv.cpp
View File

@ -35,7 +35,7 @@ PairOxdna2Excv::~PairOxdna2Excv()
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3],
void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3],
double /*e3*/[3], double rs[3], double rb[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408, d_cb=+0.4;

2
src/USER-CGDNA/pair_oxdna_coaxstk.cpp
View File

@ -59,7 +59,7 @@ PairOxdnaCoaxstk::~PairOxdnaCoaxstk()
memory->destroy(cut_cxst_hi);
memory->destroy(cut_cxst_lc);
memory->destroy(cut_cxst_hc);
memory->destroy(cutsq_cxst_hc);
memory->destroy(cutsq_cxst_hc);
memory->destroy(b_cxst_lo);
memory->destroy(b_cxst_hi);

2
src/USER-CGDNA/pair_oxdna_excv.cpp
View File

@ -86,7 +86,7 @@ PairOxdnaExcv::~PairOxdnaExcv()
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxDNA
------------------------------------------------------------------------- */
void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3],
void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3],
double /*e3*/[3], double rs[3], double rb[3])
{
double d_cs=-0.4, d_cb=+0.4;

2
src/USER-CGDNA/pair_oxdna_hbond.cpp
View File

@ -42,7 +42,7 @@ PairOxdnaHbond::PairOxdnaHbond(LAMMPS *lmp) : Pair(lmp)
// sequence-specific base-pairing strength
// A:0 C:1 G:2 T:3, 5'- [i][j] -3'
alpha_hb[0][0] = 1.00000;
alpha_hb[0][1] = 1.00000;
alpha_hb[0][2] = 1.00000;

4
src/USER-CGDNA/pair_oxdna_stk.cpp
View File

@ -43,7 +43,7 @@ PairOxdnaStk::PairOxdnaStk(LAMMPS *lmp) : Pair(lmp)
// sequence-specific stacking strength
// A:0 C:1 G:2 T:3, 5'- [i][j] -3'
eta_st[0][0] = 1.11960;
eta_st[0][0] = 1.11960;
eta_st[0][1] = 1.00852;
eta_st[0][2] = 0.96950;
eta_st[0][3] = 0.99632;
@ -121,7 +121,7 @@ PairOxdnaStk::~PairOxdnaStk()
tally energy and virial into global and per-atom accumulators
NOTE: Although this is a pair style interaction, the algorithm below
follows the virial incrementation of the bond style. This is because
follows the virial incrementation of the bond style. This is because
the bond topology is used in the main compute loop.
------------------------------------------------------------------------- */

2
src/USER-MEAMC/meam_setup_done.cpp
View File

@ -454,7 +454,7 @@ MEAM::phi_meam(double r, int a, int b)
F1 = embedding(this->A_meam[a], this->Ec_meam[a][a], rhobar1, dF);
F2 = embedding(this->A_meam[b], this->Ec_meam[b][b], rhobar2, dF);
// compute Rose function, I.16
Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],

2
src/USER-MISC/compute_gyration_shape.cpp
View File

@ -39,7 +39,7 @@ ComputeGyrationShape::ComputeGyrationShape(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 0;
// ID of compute gyration
// ID of compute gyration
int n = strlen(arg[3]) + 1;
id_gyration = new char[n];
strcpy(id_gyration,arg[3]);

58
src/USER-MISC/compute_hma.cpp
View File

@ -78,21 +78,21 @@ using namespace LAMMPS_NS;
ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), id_temp(NULL), deltaR(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute hma command");
if (narg < 4) error->all(FLERR,"Illegal compute hma command");
if (igroup) error->all(FLERR,"Compute hma must use group all");
if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}
if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
}
create_attribute = 1;
extscalar = 1;
timeflag = 1;
// (from compute displace/atom) create a new fix STORE style
// our new fix's id (id_fix)= compute-ID + COMPUTE_STORE
create_attribute = 1;
extscalar = 1;
timeflag = 1;
// (from compute displace/atom) create a new fix STORE style
// our new fix's id (id_fix)= compute-ID + COMPUTE_STORE
// our new fix's group = same as compute group
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
@ -100,30 +100,30 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[6];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
// calculate xu,yu,zu for fix store array
// skip if reset from restart file
if (fix->restart_reset) fix->restart_reset = 0;
if (fix->restart_reset) fix->restart_reset = 0;
else {
double **xoriginal = fix->astore;
double **xoriginal = fix->astore;
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
domain->unmap(x[i],image[i],xoriginal[i]);
domain->unmap(x[i],image[i],xoriginal[i]);
}
vector_flag = 1;
@ -175,7 +175,7 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
memory->create(vector, size_vector, "hma:vector");
if (computeU>-1 || computeCv>-1) {
peflag = 1;
peflag = 1;
}
if (computeP>-1) {
pressflag = 1;
@ -209,9 +209,9 @@ void ComputeHMA::init() {
}
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
}
void ComputeHMA::init_list(int /* id */, NeighList *ptr)
@ -224,22 +224,22 @@ void ComputeHMA::setup()
int dummy=0;
int ifix = modify->find_fix(id_temp);
if (ifix < 0) error->all(FLERR,"Could not find compute hma temperature ID");
double * temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);
double * temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);
if (temperat==NULL) error->all(FLERR,"Could not find compute hma temperature ID");
finaltemp = * temperat;
finaltemp = * temperat;
// set fix which stores original atom coords
int ifix2 = modify->find_fix(id_fix);
if (ifix2 < 0) error->all(FLERR,"Could not find hma store fix ID");
fix = (FixStore *) modify->fix[ifix2];
fix = (FixStore *) modify->fix[ifix2];
}
/* ---------------------------------------------------------------------- */
void ComputeHMA::compute_vector()
{
invoked_vector = update->ntimestep;
invoked_vector = update->ntimestep;
// grow deltaR array if necessary
if (comm_forward>0 && atom->nmax > nmax) {
@ -257,7 +257,7 @@ void ComputeHMA::compute_vector()
int nlocal = atom->nlocal;
double *h = domain->h;
double xprd = domain->xprd;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;

2
src/USER-MISC/compute_hma.h
View File

@ -64,4 +64,4 @@ class ComputeHMA : public Compute {
#endif
#endif

6
src/USER-MISC/pair_cosine_squared.cpp
View File

@ -125,7 +125,7 @@ void PairCosineSquared::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6)
error->all(FLERR, "Incorrect args for pair coefficients (too few or too many)");
if (!allocated)
allocate();
@ -459,7 +459,7 @@ double PairCosineSquared::single(int /* i */, int /* j */, int itype, int jtype,
double &fforce)
{
double r, r2inv, r6inv, cosone, force, energy;
r = sqrt(rsq);
if (r <= sigma[itype][jtype]) {
@ -478,7 +478,7 @@ double PairCosineSquared::single(int /* i */, int /* j */, int itype, int jtype,
}
} else {
cosone = cos(MY_PI*(r-sigma[itype][jtype]) / (2.0*w[itype][jtype]));
force = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) *
force = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) *
sin(MY_PI*(r-sigma[itype][jtype]) / w[itype][jtype]) / r;
energy = -epsilon[itype][jtype]*cosone*cosone;
}

2
src/USER-MISC/pair_extep.cpp
View File

@ -757,7 +757,7 @@ void PairExTeP::read_file(char *file)
// skip line if it is a leftover from the previous section,
// which can be identified by having 3 elements (instead of 2)
// as first words.
if (isupper(words[0][0]) && isupper(words[1][0]) && isupper(words[2][0]))
continue;

4
src/USER-MISC/pair_ilp_graphene_hbn.cpp
View File

@ -442,7 +442,7 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
van der Waals forces and energy
------------------------------------------------------------------------- */
@ -540,7 +540,7 @@ void PairILPGrapheneHBN::calc_FvdW(int eflag, int /* vflag */)
}
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */

6
src/USER-MISC/pair_kolmogorov_crespi_full.cpp
View File

@ -444,7 +444,7 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
van der Waals forces and energy
------------------------------------------------------------------------- */
@ -540,7 +540,7 @@ void PairKolmogorovCrespiFull::calc_FvdW(int eflag, int /* vflag */)
}
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
@ -790,7 +790,7 @@ void PairKolmogorovCrespiFull::calc_normal()
memory->create(dnormal,3,3,3,nmax,"KolmogorovCrespiFull:dnormal");
}
inum = list->inum;
inum = list->inum;
ilist = list->ilist;
//Calculate normals
for (ii = 0; ii < inum; ii++) {

230
src/USER-MISC/pair_local_density.cpp
View File

@ -61,9 +61,9 @@ static const char cite_pair_local_density[] =
PairLocalDensity::PairLocalDensity(LAMMPS *lmp) : Pair(lmp)
{
restartinfo = 0;
one_coeff = 1;
one_coeff = 1;
single_enable = 1;
// stuff read from tabulated file
nLD = 0;
nrho = 0;
@ -81,14 +81,14 @@ PairLocalDensity::PairLocalDensity(LAMMPS *lmp) : Pair(lmp)
lowercutsq = NULL;
frho = NULL;
rho = NULL;
// splined arrays
frho_spline = NULL;
// per-atom arrays
nmax = 0;
fp = NULL;
localrho = NULL;
localrho = NULL;
// set comm size needed by this pair
comm_forward = 1;
@ -114,10 +114,10 @@ PairLocalDensity::~PairLocalDensity()
}
memory->destroy(frho_spline);
memory->destroy(rho_min);
memory->destroy(rho_min);
memory->destroy(rho_max);
memory->destroy(delta_rho);
memory->destroy(delta_rho);
memory->destroy(c0);
memory->destroy(c2);
memory->destroy(c4);
@ -137,37 +137,37 @@ PairLocalDensity::~PairLocalDensity()
void PairLocalDensity::compute(int eflag, int vflag)
{
int i,j,ii,jj,m,k,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double rsqinv, phi, uLD, dphi, evdwl,fpair;
double p, *coeff;
int *ilist,*jlist,*numneigh,**firstneigh;
phi = uLD = evdwl = fpair = rsqinv = 0.0;
phi = uLD = evdwl = fpair = rsqinv = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
/* localrho = LD at each atom
fp = derivative of embedding energy at each atom for each LD potential
uLD = embedding energy of each atom due to each LD potential*/
uLD = embedding energy of each atom due to each LD potential*/
// grow LD and fp arrays if necessary
// need to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(localrho);
memory->destroy(fp);
nmax = atom->nmax;
nmax = atom->nmax;
memory->create(localrho, nLD, nmax, "pairLD:localrho");
memory->create(fp, nLD, nmax, "pairLD:fp");
}
double **x = atom->x;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
@ -179,13 +179,13 @@ void PairLocalDensity::compute(int eflag, int vflag)
if (newton_pair) {
m = nlocal + atom->nghost;
for (k = 0; k < nLD; k++) {
for (i = 0; i < m; i++) {
for (k = 0; k < nLD; k++) {
for (i = 0; i < m; i++) {
localrho[k][i] = 0.0;
fp[k][i] = 0.0;
}
}
}
}
}
else {
for (k = 0; k < nLD; k++){
for (i = 0; i < nlocal; i++) {
@ -196,7 +196,7 @@ void PairLocalDensity::compute(int eflag, int vflag)
}
// loop over neighs of central atoms and types of LDs
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -205,19 +205,19 @@ void PairLocalDensity::compute(int eflag, int vflag)
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
jtype = type[j];
// calculate distance-squared between i,j atom-types
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx*delx + dely*dely + delz*delz;
// calculating LDs based on central and neigh filters
for (k = 0; k < nLD; k++) {
@ -230,36 +230,36 @@ void PairLocalDensity::compute(int eflag, int vflag)
else {
phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
}
localrho[k][i] += (phi * b[k][jtype]);
/*checking for both i,j is necessary
localrho[k][i] += (phi * b[k][jtype]);
/*checking for both i,j is necessary
since a half neighbor list is processed.*/
if (newton_pair || j<nlocal) {
localrho[k][j] += (phi * b[k][itype]);
localrho[k][j] += (phi * b[k][itype]);
}
}
}
}
}
// communicate and sum LDs over all procs
if (newton_pair) comm->reverse_comm_pair(this);
//
//
for (ii = 0; ii < inum; ii++) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
uLD = 0.0;
uLD = 0.0;
for (k = 0; k < nLD; k++) {
/*skip over this loop if the LD potential
/*skip over this loop if the LD potential
is not intendend for central atomtype <itype>*/
if (!(a[k][itype])) continue;
if (!(a[k][itype])) continue;
// linear extrapolation at rho_min and rho_max
if (localrho[k][i] <= rho_min[k]) {
coeff = frho_spline[k][0];
fp[k][i] = coeff[2];
@ -284,14 +284,14 @@ void PairLocalDensity::compute(int eflag, int vflag)
if (eflag) {
if (eflag_global) eng_vdwl += uLD;
if (eflag_atom) eatom[i] += uLD;
if (eflag_atom) eatom[i] += uLD;
}
}
// communicate LD and fp to all procs
comm->forward_comm_pair(this);
// compute forces on each atom
// loop over neighbors of my atoms
@ -306,7 +306,7 @@ void PairLocalDensity::compute(int eflag, int vflag)
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
@ -316,19 +316,19 @@ void PairLocalDensity::compute(int eflag, int vflag)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
// calculate force between two atoms
// calculate force between two atoms
fpair = 0.0;
if (rsq < cutforcesq) { // global cutoff check
rsqinv = 1.0/rsq;
for (k = 0; k < nLD; k++) {
if (rsq >= lowercutsq[k] && rsq < uppercutsq[k]) {
dphi = rsq * (2.0*c2[k] + rsq * (4.0*c4[k] + 6.0*c6[k]*rsq));
fpair += -(a[k][itype]*b[k][jtype]*fp[k][i] + a[k][jtype]*b[k][itype]*fp[k][j]) * dphi;
fpair += -(a[k][itype]*b[k][jtype]*fp[k][i] + a[k][jtype]*b[k][itype]*fp[k][j]) * dphi;
}
}
fpair *= rsqinv;
}
fpair *= rsqinv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
@ -337,19 +337,19 @@ void PairLocalDensity::compute(int eflag, int vflag)
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
/*eng_vdwl has already been completely built,
/*eng_vdwl has already been completely built,
so no need to add anything here*/
if (eflag) evdwl = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
@ -362,7 +362,7 @@ void PairLocalDensity::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(setflag,n+1,n+1,"pair:setflag");
@ -430,7 +430,7 @@ void PairLocalDensity::init_style()
// request half neighbor list
array2spline();
// half neighbor request
neighbor->request(this);
}
@ -446,7 +446,7 @@ double PairLocalDensity::init_one(int /* i */, int /* j */)
cutmax = 0.0;
for (int k = 0; k < nLD; k++)
cutmax = MAX(cutmax,uppercut[k]);
cutforcesq = cutmax*cutmax;
return cutmax;
@ -454,7 +454,7 @@ double PairLocalDensity::init_one(int /* i */, int /* j */)
/*--------------------------------------------------------------------------
pair_write functionality for this pair style that gives just a snap-shot
pair_write functionality for this pair style that gives just a snap-shot
of the LD potential without doing an actual MD run
---------------------------------------------------------------------------*/
@ -473,7 +473,7 @@ double PairLocalDensity::single(int /* i */, int /* j */, int itype, int jtype,
for (k = 0; k < nLD; k++) {
LD[k][1] = 0.0; // itype:- 1
LD[k][2] = 0.0; // jtype:- 2
}
}
rsqinv = 1.0/rsq;
for (k = 0; k < nLD; k++) {
@ -487,13 +487,13 @@ double PairLocalDensity::single(int /* i */, int /* j */, int itype, int jtype,
phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
}
LD[k][1] += (phi * b[k][jtype]);
LD[k][2] += (phi * b[k][itype]);
LD[k][2] += (phi * b[k][itype]);
}
for (k = 0; k < nLD; k++) {
if (a[k][itype]) index = 1;
if (a[k][jtype]) index = 2;
if (LD[k][index] <= rho_min[k]) {
coeff = frho_spline[k][0];
dFdrho = coeff[2];
@ -545,43 +545,43 @@ void PairLocalDensity::array2spline() {
}
/* ----------------------------------------------------------------------
(one-dimensional) cubic spline interpolation sub-routine,
which determines the coeffs for a clamped cubic spline
/* ----------------------------------------------------------------------
(one-dimensional) cubic spline interpolation sub-routine,
which determines the coeffs for a clamped cubic spline
given tabulated data
------------------------------------------------------------------------*/
void PairLocalDensity::interpolate_cbspl(int n, double delta,
double *f, double **spline)
void PairLocalDensity::interpolate_cbspl(int n, double delta,
double *f, double **spline)
{
/* inputs:
n number of interpolating points
f array containing function values to
be interpolated; f[i] is the function
value corresponding to x[i]
('x' refers to the independent var)
delta difference in tabulated values of x
outputs: (packaged as columns of the coeff matrix)
coeff_b coeffs of linear terms
coeff_c coeffs of quadratic terms
coeff_d coeffs of cubic terms
spline matrix that collects b,c,d
other parameters:
fpa derivative of function at x=a
fpb derivative of function at x=b
*/
double *dl, *dd, *du;
double *coeff_b, *coeff_c, *coeff_d;
double fpa, fpb;
int i;
coeff_b = new double [n];
coeff_c = new double [n];
coeff_d = new double [n];
@ -598,11 +598,11 @@ void PairLocalDensity::interpolate_cbspl(int n, double delta,
// set slopes at beginning and end
fpa = 0.;
fpb = 0.;
for ( i = 0; i < n-1; i++ ) {
dl[i] = du[i] = delta;
}
dd[0] = 2.0 * delta;
dd[n-1] = 2.0 * delta;
coeff_c[0] = ( 3.0 / delta ) * ( f[1] - f[0] ) - 3.0 * fpa;
@ -612,20 +612,20 @@ void PairLocalDensity::interpolate_cbspl(int n, double delta,
coeff_c[i+1] = ( 3.0 / delta ) * ( f[i+2] - f[i+1] ) -
( 3.0 / delta ) * ( f[i+1] - f[i] );
}
// tridiagonal solver
for ( i = 0; i < n-1; i++ ) {
du[i] /= dd[i];
dd[i+1] -= dl[i]*du[i];
}
coeff_c[0] /= dd[0];
for ( i = 1; i < n; i++ )
coeff_c[i] = ( coeff_c[i] - dl[i-1] * coeff_c[i-1] ) / dd[i];
for ( i = n-2; i >= 0; i-- )
coeff_c[i] -= coeff_c[i+1] * du[i];
for ( i = 0; i < n-1; i++ ) {
coeff_d[i] = ( coeff_c[i+1] - coeff_c[i] ) / ( 3.0 * delta );
coeff_b[i] = ( f[i+1] - f[i] ) / delta - delta * ( coeff_c[i+1] + 2.0*coeff_c[i] ) / 3.0;
@ -648,7 +648,7 @@ void PairLocalDensity::interpolate_cbspl(int n, double delta,
spline[i][1] = 2.0*spline[i][4]/delta;
spline[i][0] = 3.0*spline[i][3]/delta;
}
delete [] coeff_b;
delete [] coeff_c;
delete [] coeff_d;
@ -662,7 +662,7 @@ void PairLocalDensity::interpolate_cbspl(int n, double delta,
------------------------------------------------------------------------- */
void PairLocalDensity::parse_file(char *filename) {
int k, n;
int me = comm->me;
FILE *fptr;
@ -680,23 +680,23 @@ void PairLocalDensity::parse_file(char *filename) {
}
double *ftmp; // tmp var to extract the complete 2D frho array from file
// broadcast number of LD potentials and number of (rho,frho) pairs
if (me == 0) {
// first 2 comment lines ignored
// first 2 comment lines ignored
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
// extract number of potentials and number of (frho, rho) points
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
sscanf(line, "%d %d", &nLD, &nrho);
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
}
MPI_Bcast(&nLD,1,MPI_INT,0,world);
MPI_Bcast(&nrho,1,MPI_INT,0,world);
// setting up all arrays to be read from files and broadcasted
memory->create(uppercut, nLD, "pairLD:uppercut");
memory->create(lowercut, nLD, "pairLD:lowercut");
@ -706,14 +706,14 @@ void PairLocalDensity::parse_file(char *filename) {
memory->create(c2, nLD, "pairLD:c2");
memory->create(c4, nLD, "pairLD:c4");
memory->create(c6, nLD, "pairLD:c6");
memory->create(rho_min, nLD, "pairLD:rho_min");
memory->create(rho_min, nLD, "pairLD:rho_min");
memory->create(rho_max, nLD, "pairLD:rho_max");
memory->create(delta_rho, nLD,"pairLD:delta_rho");
memory->create(ftmp, nrho*nLD, "pairLD:ftmp");
// setting up central and neighbor atom filters
// setting up central and neighbor atom filters
memory->create(a, nLD, atom->ntypes+1 , "pairLD:a");
memory->create(b, nLD, atom->ntypes+1, "pairLD:b");
memory->create(b, nLD, atom->ntypes+1, "pairLD:b");
if (me == 0) {
for (n = 1; n <= atom->ntypes; n++){
for (k = 0; k < nLD; k++) {
@ -721,17 +721,17 @@ void PairLocalDensity::parse_file(char *filename) {
b[k][n] = 0;
}
}
}
}
// read file block by block
if (me == 0) {
for (k = 0; k < nLD; k++) {
// parse upper and lower cut values
// parse upper and lower cut values
if (fgets(line,MAXLINE,fptr)==NULL) break;
sscanf(line, "%lf %lf", &lowercut[k], &uppercut[k]);
// parse and broadcast central atom filter
utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
char *tmp = strtok(line, " /t/n/r/f");
@ -739,27 +739,27 @@ void PairLocalDensity::parse_file(char *filename) {
a[k][atoi(tmp)] = 1;
tmp = strtok(NULL, " /t/n/r/f");
}
// parse neighbor atom filter
utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
tmp = strtok(line, " /t/n/r/f");
while (tmp != NULL) {
while (tmp != NULL) {
b[k][atoi(tmp)] = 1;
tmp = strtok(NULL, " /t/n/r/f");
}
// parse min, max and delta rho values
utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
sscanf(line, "%lf %lf %lf", &rho_min[k], &rho_max[k], &delta_rho[k]);
// recompute delta_rho from scratch for precision
delta_rho[k] = (rho_max[k] - rho_min[k]) / (nrho - 1);
// parse tabulated frho values from each line into temporary array
for (n = 0; n < nrho; n++) {
for (n = 0; n < nrho; n++) {
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
sscanf(line, "%lf", &ftmp[k*nrho + n]);
}
// ignore blank line at the end of every block
utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
@ -778,7 +778,7 @@ void PairLocalDensity::parse_file(char *filename) {
}
}
// Broadcast all parsed arrays
// Broadcast all parsed arrays
MPI_Bcast(&lowercut[0], nLD, MPI_DOUBLE, 0, world);
MPI_Bcast(&uppercut[0], nLD, MPI_DOUBLE, 0, world);
MPI_Bcast(&lowercutsq[0], nLD, MPI_DOUBLE, 0, world);
@ -800,8 +800,8 @@ void PairLocalDensity::parse_file(char *filename) {
// set up rho and frho arrays
memory->create(rho, nLD, nrho, "pairLD:rho");
memory->create(frho, nLD, nrho, "pairLD:frho");
memory->create(frho, nLD, nrho, "pairLD:frho");
for (k = 0; k < nLD; k++) {
for (n = 0; n < nrho; n++) {
rho[k][n] = rho_min[k] + n*delta_rho[k];
@ -812,7 +812,7 @@ void PairLocalDensity::parse_file(char *filename) {
// delete temporary array
memory->destroy(ftmp);
}
/* ----------------------------------------------------------------------
communication routines
------------------------------------------------------------------------- */
@ -820,16 +820,16 @@ void PairLocalDensity::parse_file(char *filename) {
int PairLocalDensity::pack_comm(int n, int *list, double *buf,
int /* pbc_flag */, int * /* pbc */) {
int i,j,k;
int m;
int m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
j = list[i];
for (k = 0; k < nLD; k++) {
buf[m++] = fp[k][j];
}
buf[m++] = fp[k][j];
}
}
return nLD;
}
@ -838,14 +838,14 @@ int PairLocalDensity::pack_comm(int n, int *list, double *buf,
void PairLocalDensity::unpack_comm(int n, int first, double *buf) {
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
for (k = 0; k < nLD; k++) {
fp[k][i] = buf[m++];
}
}
}
}
/* ---------------------------------------------------------------------- */
@ -876,7 +876,7 @@ void PairLocalDensity::unpack_reverse_comm(int n, int *list, double *buf) {
j = list[i];
for (k = 0; k < nLD; k++) {
localrho[k][j] += buf[m++];
}
}
}
}

22
src/USER-MISC/pair_local_density.h
View File

@ -13,8 +13,8 @@
pair_LocalDensity written by:
Tanmoy Sanyal and M. Scott Shell from UC Santa Barbara
David Rosenberger: TU Darmstadt
-------------------------------------------------------------------------*/
-------------------------------------------------------------------------*/
#ifdef PAIR_CLASS
@ -40,7 +40,7 @@ class PairLocalDensity : public Pair {
void init_style();
double init_one(int, int);
double single(int, int, int, int, double, double, double, double &);
virtual int pack_comm(int, int *, double *, int, int *);
virtual void unpack_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
@ -51,7 +51,7 @@ class PairLocalDensity : public Pair {
protected:
//------------------------------------------------------------------------
//This information is read from the tabulated input file
int nLD, nrho; // number of LD types
int **a, **b; // central and neigh atom filters
double *uppercut, *lowercut; // upper and lower cutoffs
@ -59,25 +59,25 @@ class PairLocalDensity : public Pair {
double *c0, *c2, *c4, *c6; // coeffs for indicator function
double *rho_min, *rho_max, *delta_rho; // min, max & grid-size for LDs
double **rho, **frho; // LD and LD function tables
//------------------------------------------------------------------------
double ***frho_spline; // splined LD potentials
double cutmax; // max cutoff for all elements
double cutforcesq; // square of global upper cutoff
int nmax; // max size of per-atom arrays
double **localrho; // per-atom LD
double **fp; // per-atom LD potential function derivative
void allocate();
// read tabulated input file
void parse_file(char *);
// convert array to spline
void array2spline();
// cubic spline interpolation
void interpolate_cbspl(int, double, double *, double **);
};

2
src/USER-MOLFILE/reader_molfile.cpp
View File

@ -282,7 +282,7 @@ bigint ReaderMolfile::read_header(double box[3][3], int &boxinfo, int &triclinic
}
// if no field info requested, just return
if (!fieldinfo) return natoms;
memory->create(fieldindex,nfield,"read_dump:fieldindex");

2
src/USER-PHONON/dynamical_matrix.cpp
View File

@ -259,7 +259,7 @@ void DynamicalMatrix::calculateMatrix()
fprintf(screen," Atoms in group = " BIGINT_FORMAT "\n", gcount);
fprintf(screen," Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
}
// emit dynlen rows of dimalpha*dynlen*dimbeta elements
update->nsteps = 0;

2
src/USER-PLUMED/fix_plumed.cpp
View File

@ -411,7 +411,7 @@ void FixPlumed::post_force(int /* vflag */)
// pass all pointers to plumed:
p->cmd("setStep",&step);
int plumedStopCondition=0;
int plumedStopCondition=0;
p->cmd("setStopFlag",&plumedStopCondition);
p->cmd("setPositions",&atom->x[0][0]);
p->cmd("setBox",&box[0][0]);

2
src/USER-REAXC/reaxc_ffield.cpp
View File

@ -154,7 +154,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Sanity checks */
if (c == 2 && !lgflag)
control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");
control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");
if (c < 9) {
snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);

6
src/USER-YAFF/pair_lj_switch3_coulgauss_long.cpp
View File

@ -133,7 +133,7 @@ void PairLJSwitch3CoulGaussLong::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
@ -202,7 +202,7 @@ void PairLJSwitch3CoulGaussLong::compute(int eflag, int vflag)
if (r>cut_lj[itype][jtype]-truncw) {
trx = (cut_lj[itype][jtype]-r)*truncwi;
tr = trx*trx*(3.0-2.0*trx);
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
forcelj = forcelj*tr + evdwl*ftr;
evdwl *= tr;
}
@ -683,7 +683,7 @@ double PairLJSwitch3CoulGaussLong::single(int i, int j, int itype, int jtype,
if (r>cut_lj[itype][jtype]-truncw) {
trx = (cut_lj[itype][jtype]-r)*truncwi;
tr = trx*trx*(3.0-2.0*trx);
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
forcelj = forcelj*tr + evdwl*ftr;
evdwl *= tr;
}

6
src/USER-YAFF/pair_mm3_switch3_coulgauss_long.cpp
View File

@ -133,7 +133,7 @@ void PairMM3Switch3CoulGaussLong::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
@ -204,7 +204,7 @@ void PairMM3Switch3CoulGaussLong::compute(int eflag, int vflag)
if (r>cut_lj[itype][jtype]-truncw) {
trx = (cut_lj[itype][jtype]-r)*truncwi;
tr = trx*trx*(3.0-2.0*trx);
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
forcelj = forcelj*tr + evdwl*ftr;
evdwl *= tr;
}
@ -683,7 +683,7 @@ double PairMM3Switch3CoulGaussLong::single(int i, int j, int itype, int jtype,
if (r>cut_lj[itype][jtype]-truncw) {
trx = (cut_lj[itype][jtype]-r)*truncwi;
tr = trx*trx*(3.0-2.0*trx);
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
ftr = 6.0*trx*(1.0-trx)*r*truncwi;
forcelj = forcelj*tr + evdwl*ftr;
evdwl *= tr;
}

2
src/comm.cpp
View File

@ -673,7 +673,7 @@ double Comm::get_comm_cutoff()
// cutoff was given and no pair style present. Otherwise print a
// warning, if the estimated bond based cutoff is larger than what
// is currently used.
if (!force->pair && (cutghostuser == 0.0)) {
maxcommcutoff = MAX(maxcommcutoff,maxbondcutoff);
} else {

2
src/compute_bond_local.cpp
View File

@ -130,7 +130,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
singleflag = 0;
velflag = 0;
for (int i = 0; i < nvalues; i++) {
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX ||
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX ||
bstyle[i] == FY || bstyle[i] == FZ) singleflag = 1;
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;

10
src/compute_orientorder_atom.cpp
View File

@ -493,7 +493,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
}
// calculate Q_l
// NOTE: optional W_l_hat and components of Q_qlcomp use these stored Q_l values
// NOTE: optional W_l_hat and components of Q_qlcomp use these stored Q_l values
int jj = 0;
for (int il = 0; il < nqlist; il++) {
@ -505,7 +505,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
qn[jj++] = qnormfac * sqrt(qm_sum);
}
// TODO:
// TODO:
// 1. [done]Need to allocate extra memory in qnarray[] for this option
// 2. [done]Need to add keyword option
// 3. [done]Need to caclulate Clebsch-Gordan/Wigner 3j coefficients
@ -673,7 +673,7 @@ void ComputeOrientOrderAtom::init_clebsch_gordan()
for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
bb2 = m2 - l;
m = aa2 + bb2 + l;
sum = 0.0;
for (int z = MAX(0, MAX(-aa2, bb2));
z <= MIN(l, MIN(l - aa2, l + bb2)); z++) {
@ -686,7 +686,7 @@ void ComputeOrientOrderAtom::init_clebsch_gordan()
factorial(aa2 + z) *
factorial(-bb2 + z));
}
cc2 = m - l;
sfaccg = sqrt(factorial(l + aa2) *
factorial(l - aa2) *
@ -695,7 +695,7 @@ void ComputeOrientOrderAtom::init_clebsch_gordan()
factorial(l + cc2) *
factorial(l - cc2) *
(2*l + 1));
sfac1 = factorial(3*l + 1);
sfac2 = factorial(l);
dcg = sqrt(sfac2*sfac2*sfac2 / sfac1);

4
src/fix_neigh_history.cpp
View File

@ -407,7 +407,7 @@ void FixNeighHistory::pre_exchange_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer)
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues,jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
@ -521,7 +521,7 @@ void FixNeighHistory::pre_exchange_no_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer)
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues, jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}

2
src/input.h
View File

@ -39,7 +39,7 @@ class Input : protected Pointers {
// substitute for variables in a string
int expand_args(int, char **, int, char **&); // expand args due to wildcard
void write_echo(const char *); // send text to active echo file pointers
protected:
char *command; // ptr to current command
int echo_screen; // 0 = no, 1 = yes

2
src/lammps.cpp
View File

@ -1013,7 +1013,7 @@ void _noopt LAMMPS::init_pkg_lists()
#undef REGION_CLASS
}
bool LAMMPS::is_installed_pkg(const char *pkg)
bool LAMMPS::is_installed_pkg(const char *pkg)
{
for (int i=0; installed_packages[i] != NULL; ++i)
if (strcmp(installed_packages[i],pkg) == 0) return true;

2
src/min.cpp
View File

@ -898,7 +898,7 @@ double Min::total_torque()
MPI_Allreduce(&ftotsqone,&ftotsqall,1,MPI_DOUBLE,MPI_SUM,world);
// multiply it by hbar so that units are in eV
return sqrt(ftotsqall) * hbar;
}

6
src/min.h
View File

@ -47,8 +47,8 @@ class Min : protected Pointers {
// methods for spin minimizers
double total_torque();
double inf_torque();
double max_torque();
double inf_torque();
double max_torque();
virtual void init_style() {}
virtual void setup_style() = 0;
@ -65,7 +65,7 @@ class Min : protected Pointers {
int external_force_clear; // clear forces locally or externally
double dmax; // max dist to move any atom in one step
int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero
int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero
// 3 = spin_cubic, 4 = spin_none
int normstyle; // TWO, MAX or INF flag for force norm evaluation

4
src/min_cg.cpp
View File

@ -91,7 +91,7 @@ int MinCG::iterate(int maxiter)
dot[0] += fvec[i]*fvec[i];
dot[1] += fvec[i]*g[i];
}
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
fatom = fextra_atom[m];
@ -119,7 +119,7 @@ int MinCG::iterate(int maxiter)
if (fmax < update->ftol*update->ftol) return FTOL;
} else if (normstyle == TWO) { // Euclidean force 2-norm
if (dotall[0] < update->ftol*update->ftol) return FTOL;
} else error->all(FLERR,"Illegal min_modify command");
} else error->all(FLERR,"Illegal min_modify command");
// update new search direction h from new f = -Grad(x) and old g
// this is Polak-Ribieri formulation

6
src/neighbor.cpp
View File

@ -1941,7 +1941,7 @@ int Neighbor::decide()
conservative shrink procedure:
compute distance each of 8 corners of box has moved since last reneighbor
reduce skin distance by sum of 2 largest of the 8 values
if reduced skin distance is negative, set to zero
if reduced skin distance is negative, set to zero
new trigger = 1/2 of reduced skin distance
for orthogonal box, only need 2 lo/hi corners
for triclinic, need all 8 corners since deformations can displace all 8
@ -1963,7 +1963,7 @@ int Neighbor::check_distance()
delz = bboxhi[2] - boxhi_hold[2];
delta2 = sqrt(delx*delx + dely*dely + delz*delz);
delta = 0.5 * (skin - (delta1+delta2));
if (delta < 0.0) delta = 0.0;
if (delta < 0.0) delta = 0.0;
deltasq = delta*delta;
} else {
domain->box_corners();
@ -1977,7 +1977,7 @@ int Neighbor::check_distance()
else if (delta > delta2) delta2 = delta;
}
delta = 0.5 * (skin - (delta1+delta2));
if (delta < 0.0) delta = 0.0;
if (delta < 0.0) delta = 0.0;
deltasq = delta*delta;
}
} else deltasq = triggersq;

2
src/read_data.cpp
View File

@ -2151,7 +2151,7 @@ void ReadData::parse_coeffs(char *line, const char *addstr,
// to avoid segfaults on empty lines
if (narg == 0) return;
if (noffset) {
int value = force->inumeric(FLERR,arg[0]);
sprintf(argoffset1,"%d",value+offset);

2
src/read_dump.cpp
View File

@ -503,7 +503,7 @@ void ReadDump::header(int fieldinfo)
yhi = box[1][1];
zlo = box[2][0];
zhi = box[2][1];
if (triclinic_snap) {
xy = box[0][2];
xz = box[1][2];

2
src/reader_native.cpp
View File

@ -113,7 +113,7 @@ void ReaderNative::skip()
only called by proc 0
------------------------------------------------------------------------- */
bigint ReaderNative::read_header(double box[3][3], int &boxinfo, int &triclinic,
bigint ReaderNative::read_header(double box[3][3], int &boxinfo, int &triclinic,
int fieldinfo, int nfield,
int *fieldtype, char **fieldlabel,
int scaleflag, int wrapflag, int &fieldflag,