Merge pull request #2813 from akohlmey/collected-small-changes

Collected small changes
This commit is contained in:
Axel Kohlmeyer 2021-06-28 20:14:01 -04:00 committed by GitHub
commit 105c86399b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 161 additions and 171 deletions

View File

@ -1,5 +1,7 @@
set(MOLFILE_INCLUDE_DIR "${LAMMPS_LIB_SOURCE_DIR}/molfile" CACHE STRING "Path to VMD molfile plugin headers")
set(MOLFILE_INCLUDE_DIRS "${MOLFILE_INCLUDE_DIR}")
add_library(molfile INTERFACE)
target_include_directories(molfile INTERFACE ${MOLFILE_INCLUDE_DIRS})
target_include_directories(molfile INTERFACE ${MOLFILE_INCLUDE_DIR})
if(NOT (CMAKE_SYSTEM_NAME STREQUAL "Windows"))
target_link_libraries(molfile INTERFACE ${CMAKE_DL_LIBS})
endif()
target_link_libraries(lammps PRIVATE molfile)

View File

@ -176,9 +176,9 @@ These fixes are not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
These fixes are not compatible with :doc:`fix shake <fix_shake>`.
Fix *temp/csld* is not compatible with :doc:`fix shake <fix_shake>`.
The fix can be used with dynamic groups as defined by the
These fixes can be used with dynamic groups as defined by the
:doc:`group <group>` command. Likewise it can be used with groups to
which atoms are added or deleted over time, e.g. a deposition
simulation. However, the conservation properties of the thermostat

View File

@ -26,7 +26,7 @@ Examples
if "${steps} > 1000" then quit
if "${myString} == a10" then quit
if "$x <= $y" then "print X is smaller = $x" else "print Y is smaller = $y"
if "(${eng} > 0.0) \|\| ($n < 1000)" then &
if "(${eng} > 0.0) || ($n < 1000)" then &
"timestep 0.005" &
elif $n<10000 &
"timestep 0.01" &

View File

@ -574,7 +574,7 @@ division and the modulo operator "%" are next; addition and
subtraction are next; the 4 relational operators "<", "<=", ">", and
">=" are next; the two remaining relational operators "==" and "!="
are next; then the logical AND operator "&&"; and finally the logical
OR operator "\|\|" and logical XOR (exclusive or) operator "\|\^" have the
OR operator "||" and logical XOR (exclusive or) operator "\|^" have the
lowest precedence. Parenthesis can be used to group one or more
portions of a formula and/or enforce a different order of evaluation
than what would occur with the default precedence.

View File

@ -221,7 +221,7 @@ void MLIAPDescriptorSO3::compute_descriptors(class MLIAPData *data)
void MLIAPDescriptorSO3::compute_forces(class MLIAPData *data)
{
int npairs = 0;
bigint npairs = 0;
for (int ii = 0; ii < data->nlistatoms; ii++) npairs += data->numneighs[ii];
so3ptr->spectrum_dxdr(data->nlistatoms, data->numneighs, data->jelems, wjelem, data->rij, nmax,

View File

@ -45,6 +45,7 @@ MLIAP_SO3::MLIAP_SO3(LAMMPS *lmp, double vrcut, int vlmax, int vnmax, double val
compute_ncoeff();
m_Nmax = (m_nmax + m_lmax + 1) * 10;
m_numYlms = (m_lmax + 1) * (m_lmax + 1);
m_ellpl1 = nullptr;
m_ellm1 = nullptr;
@ -95,6 +96,10 @@ MLIAP_SO3::MLIAP_SO3(LAMMPS *lmp, double vrcut, int vlmax, int vnmax, double val
m_clisttot_i = nullptr;
m_init_arrays = 0;
m_dfac_l1 = m_dfac_l2 = 0;
m_pfac_l1 = m_pfac_l2 = 0;
m_idxu_count = m_idxy_count = 0;
alloc_init = alloc_arrays = 0.0;
}
/* ---------------------------------------------------------------------- */
@ -171,7 +176,7 @@ void MLIAP_SO3::init()
memory->create(m_ellpl1, totali, "MLIAP_SO3:m_ellpl1");
memory->destroy(m_ellm1);
memory->create(m_ellm1, totali, "MLIAP_SO3:m_ellm1");
alloc_init = 2.0*totali*sizeof(double);
alloc_init = 2.0 * totali * sizeof(double);
for (int l = 1; l < m_lmax + 1; l++) {
m_ellpl1[l] = get_sum(0, l + 2, 1, 2);
@ -184,13 +189,13 @@ void MLIAP_SO3::init()
totali = m_pfac_l1 * m_pfac_l2;
memory->destroy(m_pfac);
memory->create(m_pfac, totali, "MLIAP_SO3:m_pfac");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (int l = 0; l < m_lmax + 2; l++)
for (int m = -l; m < l + 1; m++)
m_pfac[l * m_pfac_l2 + m] = sqrt((2.0 * l + 1.0) * pfac1) * powsign(m);
m_dfac_l1 = m_lmax + 1;
m_dfac_l2 = (m_lmax + 1) * (m_lmax + 1) + 1;
m_dfac_l2 = m_numYlms + 1;
totali = m_dfac_l1 * m_dfac_l2;
memory->destroy(m_dfac0);
memory->create(m_dfac0, totali, "MLIAP_SO3:m_dfac0");
@ -204,7 +209,7 @@ void MLIAP_SO3::init()
memory->create(m_dfac4, totali, "MLIAP_SO3:m_dfac4");
memory->destroy(m_dfac5);
memory->create(m_dfac5, totali, "MLIAP_SO3:m_dfac5");
alloc_init += 6.0*totali*sizeof(double);
alloc_init += 6.0 * totali * sizeof(double);
for (int l = 1; l < m_lmax + 1; l++)
for (int m = -l; m < l + 1; m++) {
@ -225,7 +230,7 @@ void MLIAP_SO3::init()
totali = m_nmax * m_nmax;
memory->destroy(m_w);
memory->create(m_w, totali, "MLIAP_SO3:w");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (i = 0; i < totali; i++) m_w[i] = 0.0;
@ -234,7 +239,7 @@ void MLIAP_SO3::init()
totali = m_nmax * m_Nmax;
memory->destroy(m_g_array);
memory->create(m_g_array, totali, "MLIAP_SO3:g_array");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (i = 0; i < totali; i++) m_g_array[i] = 0.0;
@ -242,11 +247,11 @@ void MLIAP_SO3::init()
int twolmax;
twolmax = 2 * (m_lmax + 1);
m_ldim = twolmax + 1;
int m_ldim = twolmax + 1;
totali = m_ldim * m_ldim;
memory->destroy(m_rootpq);
memory->create(m_rootpq, totali, "MLIAP_SO3:rootpq");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (i = 0; i < totali; i++) m_rootpq[i] = 0.0;
@ -255,18 +260,18 @@ void MLIAP_SO3::init()
memory->destroy(m_idxu_block);
memory->create(m_idxu_block, m_ldim, "MLIAP_SO3:idxu_bloc");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (i = 0; i < m_ldim; i++) m_idxu_block[i] = 0;
totali = square(m_lmax + 2);
memory->destroy(m_idxylm);
memory->create(m_idxylm, totali, "MLIAP_SO3:idxylm");
alloc_init += totali*sizeof(double);
alloc_init += totali * sizeof(double);
for (i = 0; i < totali; i++) m_idxylm[i] = 0;
m_idxu_count = 0, m_idxy_count = 0;
m_idxu_count = m_idxy_count = 0;
for (int l = 0; l < m_ldim; l++) {
m_idxu_block[l] = m_idxu_count;
@ -279,8 +284,6 @@ void MLIAP_SO3::init()
m_idxu_count += 1;
}
}
m_numYlms = (m_lmax + 1) * (m_lmax + 1);
}
/* ---------------------------------------------------------------------- */
@ -292,54 +295,54 @@ void MLIAP_SO3::init_arrays(int nlocal, int ncoefs)
memory->create(m_plist_r, totali, "MLIAP_SO3:m_plist_r");
memory->destroy(m_plist_i);
memory->create(m_plist_i, totali, "MLIAP_SO3:m_plist_i");
alloc_arrays = 2.0*totali*sizeof(double);
alloc_arrays = 2.0 * totali * sizeof(double);
totali = m_nmax * m_numYlms;
memory->destroy(m_clist_r);
memory->create(m_clist_r, totali, "MLIAP_SO3:m_clist_r");
memory->destroy(m_clist_i);
memory->create(m_clist_i, totali, "MLIAP_SO3:m_clist_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = m_idxu_count;
memory->destroy(m_ulist_r);
memory->create(m_ulist_r, totali, "MLIAP_SO3:m_ulist_r");
memory->destroy(m_ulist_i);
memory->create(m_ulist_i, totali, "MLIAP_SO3:m_ulist_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = (m_lmax + 2) * (m_lmax + 2);
memory->destroy(m_Ylms_r);
memory->create(m_Ylms_r, totali, "MLIAP_SO3:m_Ylms_r");
memory->destroy(m_Ylms_i);
memory->create(m_Ylms_i, totali, "MLIAP_SO3:m_Ylms_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = (m_lmax + 1) * (m_lmax + 1) * 3;
totali = m_numYlms * 3;
memory->destroy(m_dYlm_r);
memory->create(m_dYlm_r, totali, "MLIAP_SO3:m_dYlm_r");
memory->destroy(m_dYlm_i);
memory->create(m_dYlm_i, totali, "MLIAP_SO3:m_dYlm_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = m_nmax * m_numYlms * 3;
memory->destroy(m_dclist_r);
memory->create(m_dclist_r, totali, "MLIAP_SO3:m_dclist_r");
memory->destroy(m_dclist_i);
memory->create(m_dclist_i, totali, "MLIAP_SO3:m_dclist_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = 3 * m_nmax * (m_nmax + 1) * (m_lmax + 1) / 2;
memory->destroy(m_tempdp_r);
memory->create(m_tempdp_r, totali, "MLIAP_SO3:m_tempdp_r");
alloc_arrays += totali*sizeof(double);
alloc_arrays += totali * sizeof(double);
totali = m_nmax * m_numYlms;
memory->destroy(m_clisttot_r);
memory->create(m_clisttot_r, totali, "MLIAP_SO3:m_clisttot_r");
memory->destroy(m_clisttot_i);
memory->create(m_clisttot_i, totali, "MLIAP_SO3:m_clisttot_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
m_init_arrays = 1;
}
@ -374,8 +377,8 @@ void MLIAP_SO3::compute_W(int nmax, double *arr)
double *sqrtD = new double[n * n];
double *tempM = new double[n * n];
double **temparr = new double*[n];
double **tempvl = new double*[n];
double **temparr = new double *[n];
double **tempvl = new double *[n];
double *tempout = new double[n];
int info;
@ -819,54 +822,54 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem
bigint totaln = 0;
bigint totali;
double Ylm_r, Ylm_i;
int i, ti;
int weight, neighbor;
double x, y, z, r;
double r_int;
int twolmax = 2 * (lmax + 1);
int findex, gindex;
int ipair = 0;
int findex;
bigint gindex;
bigint ipair = 0;
double sfac;
findex = m_nmax * (m_lmax + 1);
for (i = 0; i < nlocal; i++) totaln += numneighs[i];
for (int i = 0; i < nlocal; i++) totaln += numneighs[i];
totali = totaln * m_Nmax * (m_lmax + 1);
memory->destroy(m_sbes_array);
memory->create(m_sbes_array, totali, "MLIAP_SO3:m_sbes_array");
memory->destroy(m_sbes_darray);
memory->create(m_sbes_darray, totali, "MLIAP_SO3:m_sbes_darray");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = totaln * m_nmax * (m_lmax + 1);
memory->destroy(m_rip_array);
memory->create(m_rip_array, totali, "MLIAP_SO3:m_rip_array");
memory->destroy(m_rip_darray);
memory->create(m_rip_darray, totali, "MLIAP_SO3:m_rip_darray");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
totali = totaln * ncoefs * 3;
memory->destroy(m_dplist_r);
memory->create(m_dplist_r, totali, "MLIAP_SO3:m_dplist_r");
memory->destroy(m_dplist_i);
memory->create(m_dplist_i, totali, "MLIAP_SO3:m_dplist_i");
alloc_arrays += 2.0*totali*sizeof(double);
alloc_arrays += 2.0 * totali * sizeof(double);
get_sbes_array(nlocal, numneighs, rij, lmax, rcut, alpha);
get_rip_array(nlocal, numneighs, rij, nmax, lmax, alpha);
totali = nlocal * ncoefs;
for (i = 0; i < totali; i++) {
totali = (bigint) nlocal * ncoefs;
for (int i = 0; i < totali; i++) {
m_plist_r[i] = 0.0;
m_plist_i[i] = 0.0;
}
for (int ii = 0; ii < nlocal; ii++) {
totali = nmax * m_numYlms;
totali = (bigint) nmax * m_numYlms;
for (ti = 0; ti < totali; ti++) {
for (bigint ti = 0; ti < totali; ti++) {
m_clisttot_r[ti] = 0.0;
m_clisttot_i[ti] = 0.0;
}
@ -883,12 +886,12 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem
r = sqrt(x * x + y * y + z * z);
if (r < SMALL) continue;
totali = nmax * m_numYlms;
for (ti = 0; ti < totali; ti++) {
totali = (bigint) nmax * m_numYlms;
for (bigint ti = 0; ti < totali; ti++) {
m_clist_r[ti] = 0.0;
m_clist_i[ti] = 0.0;
}
for (ti = 0; ti < m_idxu_count; ti++) {
for (bigint ti = 0; ti < m_idxu_count; ti++) {
m_ulist_r[ti] = 0.0;
m_ulist_i[ti] = 0.0;
}
@ -914,18 +917,17 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem
}
}
totali = nmax * m_numYlms;
for (int tn = 0; tn < totali; tn++) {
totali = (bigint) nmax * m_numYlms;
for (bigint tn = 0; tn < totali; tn++) {
m_clist_r[tn] = m_clist_r[tn] * double(weight);
m_clist_i[tn] = m_clist_i[tn] * double(weight);
}
for (int tn = 0; tn < totali; tn++) {
for (bigint tn = 0; tn < totali; tn++) {
m_clisttot_r[tn] += m_clist_r[tn];
m_clisttot_i[tn] += m_clist_i[tn];
}
}
compute_pi(nmax, lmax, m_clisttot_r, m_clisttot_i, m_numYlms, m_plist_r, m_plist_i, ncoefs, ii);
}
@ -935,7 +937,8 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem
/* ---------------------------------------------------------------------- */
void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *wjelem, double **rij,
int nmax, int lmax, double rcut, double alpha, int npairs, int ncoefs)
int nmax, int lmax, double rcut, double alpha, bigint npairs,
int ncoefs)
{
bigint totaln = 0;
bigint totali;
@ -946,15 +949,14 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
double dexpfac[3];
double dfact[6];
int i, ti;
double xcov0_r, xcov0_i, xcovpl1_r, xcovpl1_i, xcovm1_r, xcovm1_i;
double comj_i;
double r_int;
double r_int_temp;
double oneofr;
int findex, gindex;
int findex;
bigint gindex;
int numps, weight, neighbor;
@ -966,7 +968,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
findex = m_nmax * (m_lmax + 1);
for (i = 0; i < nlocal; i++) totaln += numneighs[i];
for (int i = 0; i < nlocal; i++) totaln += numneighs[i];
totali = totaln * m_Nmax * (m_lmax + 1);
memory->destroy(m_sbes_array);
@ -988,7 +990,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
totali = npairs * ncoefs * 3;
for (i = 0; i < totali; i++) {
for (int i = 0; i < totali; i++) {
m_dplist_r[i] = 0.0;
m_dplist_i[i] = 0.0;
}
@ -1004,7 +1006,7 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
for (int ii = 0; ii < nlocal; ii++) {
totali = nmax * m_numYlms;
for (ti = 0; ti < totali; ti++) {
for (bigint ti = 0; ti < totali; ti++) {
m_clisttot_r[ti] = 0.0;
m_clisttot_i[ti] = 0.0;
}
@ -1024,12 +1026,12 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
if (r < SMALL) continue;
totali = nmax * m_numYlms;
for (ti = 0; ti < totali; ti++) {
for (bigint ti = 0; ti < totali; ti++) {
m_clist_r[ti] = 0.0;
m_clist_i[ti] = 0.0;
}
for (ti = 0; ti < m_idxu_count; ti++) {
for (bigint ti = 0; ti < m_idxu_count; ti++) {
m_ulist_r[ti] = 0.0;
m_ulist_i[ti] = 0.0;
}
@ -1056,12 +1058,12 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
}
totali = nmax * m_numYlms;
for (int tn = 0; tn < totali; tn++) {
for (bigint tn = 0; tn < totali; tn++) {
m_clist_r[tn] = m_clist_r[tn] * double(weight);
m_clist_i[tn] = m_clist_i[tn] * double(weight);
}
for (int tn = 0; tn < totali; tn++) {
for (bigint tn = 0; tn < totali; tn++) {
m_clisttot_r[tn] += m_clist_r[tn];
m_clisttot_i[tn] += m_clist_i[tn];
}
@ -1080,13 +1082,13 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
r = sqrt(x * x + y * y + z * z);
if (r < SMALL) continue;
totali = nmax * m_numYlms * 3;
for (int tn = 0; tn < totali; tn++) {
totali = (bigint) nmax * m_numYlms * 3;
for (bigint tn = 0; tn < totali; tn++) {
m_dclist_r[tn] = 0.0;
m_dclist_i[tn] = 0.0;
}
for (ti = 0; ti < m_idxu_count; ti++) {
for (int ti = 0; ti < m_idxu_count; ti++) {
m_ulist_r[ti] = 0.0;
m_ulist_i[ti] = 0.0;
}
@ -1098,8 +1100,8 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
rvec[0] = x;
rvec[1] = y;
rvec[2] = z;
totali = (lmax + 2) * (lmax + 2);
for (int tn = 0; tn < totali; tn++) {
totali = ((bigint) lmax + 2) * (lmax + 2);
for (bigint tn = 0; tn < totali; tn++) {
m_Ylms_r[tn] = 0.0;
m_Ylms_i[tn] = 0.0;
}
@ -1113,8 +1115,8 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
}
}
totali = (lmax + 1) * (lmax + 1) * 3;
for (int tn = 0; tn < totali; tn++) {
totali = ((bigint) lmax + 1) * (lmax + 1) * 3;
for (bigint tn = 0; tn < totali; tn++) {
m_dYlm_r[tn] = 0.0;
m_dYlm_i[tn] = 0.0;
}
@ -1221,19 +1223,19 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w
}
/////// end compute_carray_wD //////////////////
totali = nmax * m_numYlms * 3;
for (int tn = 0; tn < totali; tn++) {
totali = (bigint) nmax * m_numYlms * 3;
for (bigint tn = 0; tn < totali; tn++) {
m_dclist_r[tn] = m_dclist_r[tn] * double(weight);
m_dclist_i[tn] = m_dclist_i[tn] * double(weight);
}
totali = numps * 3;
for (ti = 0; ti < totali; ti++) m_tempdp_r[ti] = 0.0;
totali = (bigint) numps * 3;
for (bigint ti = 0; ti < totali; ti++) m_tempdp_r[ti] = 0.0;
compute_dpidrj(nmax, lmax, m_clisttot_r, m_clisttot_i, m_numYlms, m_dclist_r, m_dclist_i,
m_numYlms, 3, m_tempdp_r, 3);
for (int tn = 0; tn < totali; tn++)
for (bigint tn = 0; tn < totali; tn++)
m_dplist_r[((idpair - 1) * (numps * 3)) + tn] += m_tempdp_r[tn];
} //for(neighbor=0;neighbor<numneighs[ii];neighbor++){

View File

@ -29,7 +29,7 @@ class MLIAP_SO3 : protected Pointers {
double *m_dfac0, *m_dfac1, *m_dfac2, *m_dfac3, *m_dfac4, *m_dfac5;
int m_dfac_l1, m_dfac_l2;
double m_rcut, m_alpha;
int m_lmax, m_nmax, m_Nmax, m_ldim;
int m_lmax, m_nmax, m_Nmax;
double *m_g_array, *m_w, *m_rootpq;
int *m_idxu_block, *m_idxylm;
int m_idxu_count, m_idxy_count;
@ -52,7 +52,7 @@ class MLIAP_SO3 : protected Pointers {
void spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem, double **rij, int nmax,
int lmax, double rcut, double alpha, int ncoefs);
void spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *wjelem, double **rij,
int nmax, int lmax, double rcut, double alpha, int npairs, int ncoefs);
int nmax, int lmax, double rcut, double alpha, bigint npairs, int ncoefs);
private:
double Cosine(double Rij, double Rc);

View File

@ -25,8 +25,8 @@ typedef Jacobi<double, double *, double **, double const *const *> Jacobi_v2;
int SO3Math::jacobin(int n, double const *const *mat, double *eval, double **evec)
{
int *midx = new int[n];
double **M = new double*[n];
double **mat_cpy = new double*[n];
double **M = new double *[n];
double **mat_cpy = new double *[n];
for (int i = 0; i < n; i++) {
mat_cpy[i] = new double[n];
@ -63,22 +63,26 @@ int SO3Math::invert_matrix(int n, double *A, double *Ainv)
for (i = 0; i < n * n; i++) Atemp[i] = A[i];
if (LUPdecompose(n, dtol, Atemp, P) != 0) return 1;
int rv = 0;
if (LUPdecompose(n, dtol, Atemp, P) == 0) {
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) b[j] = 0.0;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) b[j] = 0.0;
b[i] = 1.0;
LUPSolve(n, Atemp, b, P);
b[i] = 1.0;
LUPSolve(n, Atemp, b, P);
for (j = 0; j < n; j++) Ainv[j * n + i] = b[j];
for (j = 0; j < n; j++) Ainv[j * n + i] = b[j];
}
} else {
rv = 1;
}
delete[] P;
delete[] b;
delete[] Atemp;
return 0;
return rv;
}
int SO3Math::LUPdecompose(int n, double dtol, double *A, int *P)
@ -96,7 +100,10 @@ int SO3Math::LUPdecompose(int n, double dtol, double *A, int *P)
Atemp = fabs(A[i * n + j]);
if (Atemp > maxA) maxA = Atemp;
}
if (maxA < dtol) return 1;
if (maxA < dtol) {
delete[] normi;
return 1;
}
normi[i] = 1.0 / maxA;
}

View File

@ -38,7 +38,7 @@ using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_compute_saed_c[] =
"compute_saed command:\n\n"
"compute_saed command: doi:10.1088/0965-0393/21/5/055020\n\n"
"@Article{Coleman13,\n"
" author = {S. P. Coleman, D. E. Spearot, L. Capolungo},\n"
" title = {Virtual diffraction analysis of Ni [010] symmetric tilt grain boundaries},\n"
@ -74,7 +74,7 @@ ComputeSAED::ComputeSAED(LAMMPS *lmp, int narg, char **arg) :
extvector = 0;
// Store radiation wavelength
lambda = atof(arg[3]);
lambda = utils::numeric(FLERR,arg[3],false,lmp);
if (lambda < 0)
error->all(FLERR,"Compute SAED: Wavelength must be greater than zero");
@ -109,30 +109,30 @@ ComputeSAED::ComputeSAED(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg],"Kmax") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal Compute SAED Command");
Kmax = atof(arg[iarg+1]);
Kmax = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (Kmax / 2 < 0 || Kmax / 2 > 6)
error->all(FLERR,"Compute SAED: |K|max/2 must be between 0 and 6 ");
iarg += 2;
} else if (strcmp(arg[iarg],"Zone") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal Compute SAED Command");
Zone[0] = atof(arg[iarg+1]);
Zone[1] = atof(arg[iarg+2]);
Zone[2] = atof(arg[iarg+3]);
Zone[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
Zone[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
Zone[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
iarg += 4;
} else if (strcmp(arg[iarg],"c") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal Compute SAED Command");
c[0] = atof(arg[iarg+1]);
c[1] = atof(arg[iarg+2]);
c[2] = atof(arg[iarg+3]);
c[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
c[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
c[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (c[0] < 0 || c[1] < 0 || c[2] < 0)
error->all(FLERR,"Compute SAED: dKs must be greater than 0");
iarg += 4;
} else if (strcmp(arg[iarg],"dR_Ewald") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal Compute SAED Command");
dR_Ewald = atof(arg[iarg+1]);
dR_Ewald = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (dR_Ewald < 0)
error->all(FLERR,"Compute SAED: dR_Ewald slice must be greater than 0");
iarg += 2;
@ -170,6 +170,7 @@ ComputeSAED::ComputeSAED(LAMMPS *lmp, int narg, char **arg) :
double *prd;
double ave_inv = 0.0;
prd = domain->prd;
if (periodicity[0]) {
prd_inv[0] = 1 / prd[0];
ave_inv += prd_inv[0];
@ -203,18 +204,18 @@ ComputeSAED::ComputeSAED(LAMMPS *lmp, int narg, char **arg) :
}
}
// Find reprical spacing and integer dimensions
// Find reciprocal spacing and integer dimensions
for (int i=0; i<3; i++) {
dK[i] = prd_inv[i]*c[i];
Knmax[i] = ceil(Kmax / dK[i]);
}
// Finding the intersection of the reciprical space and Ewald sphere
// Finding the intersection of the reciprocal space and Ewald sphere
int n = 0;
double dinv2, r2, EmdR2, EpdR2;
double K[3];
// Zone flag to capture entire recrocal space volume
// Zone flag to capture entire reciprocal space volume
if ((Zone[0] == 0) && (Zone[1] == 0) && (Zone[2] == 0)) {
for (int k = -Knmax[2]; k <= Knmax[2]; k++) {
for (int j = -Knmax[1]; j <= Knmax[1]; j++) {
@ -250,13 +251,10 @@ ComputeSAED::ComputeSAED(LAMMPS *lmp, int narg, char **arg) :
}
}
if (me == 0) {
if (screen && echo) {
fprintf(screen,"-----\nCompute SAED id:%s, # of atoms:%d, # of relp:%d\n",id,natoms,n);
fprintf(screen,"Reciprocal point spacing in k1,k2,k3 = %g %g %g\n-----\n",
dK[0], dK[1], dK[2]);
}
}
if (me == 0 && echo)
utils::logmesg(lmp,"-----\nCompute SAED id:{}, # of atoms:{}, # of relp:{}\n"
"Reciprocal point spacing in k1,k2,k3 = {:.8} {:.8} {:.8}\n-----\n",
id,natoms,n,dK[0],dK[1],dK[2]);
nRows = n;
size_vector = n;
@ -347,10 +345,8 @@ void ComputeSAED::compute_vector()
{
invoked_vector = update->ntimestep;
if (me == 0 && echo) {
if (screen)
fprintf(screen,"-----\nComputing SAED intensities");
}
if (me == 0 && echo)
utils::logmesg(lmp,"-----\nComputing SAED intensities");
double t0 = MPI_Wtime();
double *Fvec = new double[2*nRows]; // Strct factor (real & imaginary)
@ -406,16 +402,10 @@ void ComputeSAED::compute_vector()
// Setting up OMP
#if defined(_OPENMP)
if (me == 0 && echo) {
if (screen)
fprintf(screen," using %d OMP threads\n",comm->nthreads);
}
if (me == 0 && echo) utils::logmesg(lmp," using {}OMP threads\n",comm->nthreads);
#endif
if (me == 0 && echo) {
if (screen)
fprintf(screen,"\n");
}
if (me == 0 && echo) utils::logmesg(lmp,"\n");
int m = 0;
double frac = 0.1;
@ -482,7 +472,7 @@ void ComputeSAED::compute_vector()
#endif
{
if (m == round(frac * nRows)) {
if (me == 0 && screen) fprintf(screen," %0.0f%% -",frac*100);
if (me == 0) utils::logmesg(lmp," {:2.0f}% -",frac*100);
frac += 0.1;
}
m++;
@ -506,10 +496,9 @@ void ComputeSAED::compute_vector()
// compute memory usage per processor
double bytes = memory_usage();
if (me == 0 && echo) {
if (screen)
fprintf(screen," 100%% \nTime elapsed during compute_saed = %0.2f sec using %0.2f Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
}
if (me == 0 && echo)
utils::logmesg(lmp," 100% \nTime elapsed during compute_saed = {:.2f} sec "
"using {:.2f} Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
delete [] xlocal;
delete [] typelocal;

View File

@ -39,7 +39,7 @@ using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_compute_xrd_c[] =
"compute_xrd command:\n\n"
"compute_xrd command: doi:10.1088/0965-0393/21/5/055020\n\n"
"@Article{Coleman13,\n"
" author = {S. P. Coleman, D. E. Spearot, L. Capolungo},\n"
" title = {Virtual diffraction analysis of Ni [010] symmetric tilt grain boundaries},\n"
@ -75,7 +75,7 @@ ComputeXRD::ComputeXRD(LAMMPS *lmp, int narg, char **arg) :
extarray = 0;
// Store radiation wavelength
lambda = atof(arg[3]);
lambda = utils::numeric(FLERR,arg[3],false,lmp);
if (lambda < 0)
error->all(FLERR,"Compute SAED: Wavelength must be greater than zero");
@ -109,33 +109,22 @@ ComputeXRD::ComputeXRD(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"2Theta") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal Compute XRD Command");
Min2Theta = atof(arg[iarg+1]) / 2;
Max2Theta = atof(arg[iarg+2]) / 2;
if (Max2Theta > MY_PI) {
Min2Theta = Min2Theta * MY_PI / 180; // converting to radians if necessary
Max2Theta = Max2Theta * MY_PI / 180;
radflag = 0;
}
if (Min2Theta <= 0)
error->all(FLERR,"Minimum 2theta value must be greater than zero");
if (Max2Theta >= MY_PI )
error->all(FLERR,"Maximum 2theta value must be less than 180 degrees");
if (Max2Theta-Min2Theta <= 0)
error->all(FLERR,"Two-theta range must be greater than zero");
Min2Theta = utils::numeric(FLERR,arg[iarg+1],false,lmp);
Max2Theta = utils::numeric(FLERR,arg[iarg+2],false,lmp);
iarg += 3;
} else if (strcmp(arg[iarg],"c") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal Compute XRD Command");
c[0] = atof(arg[iarg+1]);
c[1] = atof(arg[iarg+2]);
c[2] = atof(arg[iarg+3]);
c[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
c[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
c[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (c[0] < 0 || c[1] < 0 || c[2] < 0)
error->all(FLERR,"Compute XRD: c's must be greater than 0");
iarg += 4;
} else if (strcmp(arg[iarg],"LP") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal Compute XRD Command");
LP = atof(arg[iarg+1]);
LP = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (!(LP == 1 || LP == 0))
error->all(FLERR,"Compute XRD: LP must have value of 0 or 1");
@ -152,6 +141,20 @@ ComputeXRD::ComputeXRD(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal Compute XRD Command");
}
// error check and process min/max 2Theta values
Min2Theta /= 2.0;
Max2Theta /= 2.0;
if (Max2Theta > MY_PI) {
Min2Theta = Min2Theta * MY_PI / 180; // converting to radians if necessary
Max2Theta = Max2Theta * MY_PI / 180;
radflag = 0;
}
if (Min2Theta <= 0)
error->all(FLERR,"Minimum 2Theta value must be greater than zero");
if (Max2Theta >= MY_PI )
error->all(FLERR,"Maximum 2Theta value must be less than 180 degrees");
if (Max2Theta-Min2Theta <= 0)
error->all(FLERR,"2Theta range must be greater than zero");
Kmax = 2 * sin(Max2Theta) / lambda;
// Calculating spacing between reciprocal lattice points
@ -231,13 +234,10 @@ ComputeXRD::ComputeXRD(LAMMPS *lmp, int narg, char **arg) :
size_array_rows = nRows;
size_array_cols = 2;
if (me == 0) {
if (screen && echo) {
fprintf(screen,"-----\nCompute XRD id:%s, # of atoms:%d, # of relp:%d\n",id,natoms,nRows);
fprintf(screen,"Reciprocal point spacing in k1,k2,k3 = %g %g %g\n-----\n",
dK[0], dK[1], dK[2]);
}
}
if (me == 0 && echo)
utils::logmesg(lmp,"-----\nCompute XRD id:{}, # of atoms:{}, # of relp:{}\n"
"Reciprocal point spacing in k1,k2,k3 = {:.8} {:.8} {:.8}\n-----\n",
id,natoms,nRows,dK[0],dK[1],dK[2]);
memory->create(array,size_array_rows,size_array_cols,"xrd:array");
memory->create(store_tmp,3*size_array_rows,"xrd:store_tmp");
@ -263,9 +263,7 @@ void ComputeXRD::init()
double ang = 0.0;
double convf = 360 / MY_PI;
if (radflag ==1) {
convf = 1;
}
if (radflag ==1) convf = 1;
int n = 0;
for (int m = 0; m < mmax; m++) {
@ -300,10 +298,7 @@ void ComputeXRD::compute_array()
{
invoked_array = update->ntimestep;
if (me == 0 && echo) {
if (screen)
fprintf(screen,"-----\nComputing XRD intensities");
}
if (me == 0 && echo) utils::logmesg(lmp, "-----\nComputing XRD intensities");
double t0 = MPI_Wtime();
@ -339,19 +334,16 @@ void ComputeXRD::compute_array()
// Setting up OMP
#if defined(_OPENMP)
if (me == 0 && echo) {
if (screen)
fprintf(screen," using %d OMP threads\n",comm->nthreads);
}
if ((me == 0) && echo) utils::logmesg(lmp," using {} OMP threads\n",comm->nthreads);
#endif
if (me == 0 && echo) {
if (screen) {
fprintf(screen,"\n");
if (LP == 1)
fprintf(screen,"Applying Lorentz-Polarization Factor During XRD Calculation 2\n");
}
if ((me == 0) && echo) {
utils::logmesg(lmp,"\n");
if (LP == 1)
utils::logmesg(lmp,"Applying Lorentz-Polarization Factor During XRD Calculation 2\n");
}
int m = 0;
double frac = 0.1;
@ -430,7 +422,7 @@ void ComputeXRD::compute_array()
#endif
{
if (m == round(frac * size_array_rows)) {
if (me == 0 && screen) fprintf(screen," %0.0f%% -",frac*100);
if (me == 0) utils::logmesg(lmp," {:2.0f}% -",frac*100);
frac += 0.1;
}
m++;
@ -484,7 +476,7 @@ void ComputeXRD::compute_array()
#endif
{
if (m == round(frac * size_array_rows)) {
if (me == 0 && screen) fprintf(screen," %0.0f%% -",frac*100 );
if (me == 0) utils::logmesg(lmp," {:2.0f}% -",frac*100);
frac += 0.1;
}
m++;
@ -509,10 +501,9 @@ void ComputeXRD::compute_array()
// compute memory usage per processor
double bytes = memory_usage();
if (me == 0 && echo) {
if (screen)
fprintf(screen," 100%% \nTime elapsed during compute_xrd = %0.2f sec using %0.2f Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
}
if (me == 0 && echo)
utils::logmesg(lmp," 100% \nTime elapsed during compute_xrd = {:.2f} sec "
"using {:.2f} Mbytes/processor\n-----\n", t2-t0, bytes/1024.0/1024.0);
delete [] scratch;
delete [] Fvec;
@ -535,4 +526,3 @@ double ComputeXRD::memory_usage()
return bytes;
}