ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #1648 from athomps/add_w_l_orientorder

Browse Source 
Add w l orientorder
This commit is contained in:
Axel Kohlmeyer
2019-09-11 14:08:33 -04:00
committed by GitHub
parent df63a53788 3ec9317d03
commit ce7c3e6864
6 changed files with 729 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
doc/src/compute_orientorder_atom.txt
View File

@ -19,6 +19,8 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
{cutoff} value = distance cutoff {cutoff} value = distance cutoff
{nnn} value = number of nearest neighbors {nnn} value = number of nearest neighbors
{degrees} values = nlvalues, l1, l2,... {degrees} values = nlvalues, l1, l2,...
{wl} value = yes or no
{wl/hat} value = yes or no
{components} value = ldegree :pre {components} value = ldegree :pre
:ule :ule
@ -27,7 +29,8 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
compute 1 all orientorder/atom compute 1 all orientorder/atom
compute 1 all orientorder/atom degrees 5 4 6 8 10 12 nnn NULL cutoff 1.5 compute 1 all orientorder/atom degrees 5 4 6 8 10 12 nnn NULL cutoff 1.5
compute 1 all orientorder/atom degrees 4 6 components 6 nnn NULL cutoff 3.0 :pre compute 1 all orientorder/atom wl/hat yes
compute 1 all orientorder/atom components 6 :pre
[Description:] [Description:]
@ -48,7 +51,7 @@ neighbors of the central atom.
The angles theta and phi are the standard spherical polar angles The angles theta and phi are the standard spherical polar angles
defining the direction of the bond vector {rij}. defining the direction of the bond vector {rij}.
The second equation defines {Ql}, which is a The second equation defines {Ql}, which is a
rotationally invariant scalar quantity obtained by summing rotationally invariant non-negative amplitude obtained by summing
over all the components of degree {l}. over all the components of degree {l}.
The optional keyword {cutoff} defines the distance cutoff The optional keyword {cutoff} defines the distance cutoff
@ -63,7 +66,7 @@ specified distance cutoff are used.
The optional keyword {degrees} defines the list of order parameters to The optional keyword {degrees} defines the list of order parameters to
be computed. The first argument {nlvalues} is the number of order be computed. The first argument {nlvalues} is the number of order
parameters. This is followed by that number of integers giving the parameters. This is followed by that number of non-negative integers giving the
degree of each order parameter. Because {Q}2 and all odd-degree order degree of each order parameter. Because {Q}2 and all odd-degree order
parameters are zero for atoms in cubic crystals (see parameters are zero for atoms in cubic crystals (see
"Steinhardt"_#Steinhardt), the default order parameters are {Q}4, "Steinhardt"_#Steinhardt), the default order parameters are {Q}4,
@ -71,7 +74,20 @@ parameters are zero for atoms in cubic crystals (see
= sqrt(7/3)/8 = 0.19094.... The numerical values of all order = sqrt(7/3)/8 = 0.19094.... The numerical values of all order
parameters up to {Q}12 for a range of commonly encountered parameters up to {Q}12 for a range of commonly encountered
high-symmetry structures are given in Table I of "Mickel et high-symmetry structures are given in Table I of "Mickel et
al."_#Mickel. al."_#Mickel, and these can be reproduced with this compute
The optional keyword {wl} will output the third-order invariants {Wl}
(see Eq. 1.4 in "Steinhardt"_#Steinhardt) for the same degrees as
for the {Ql} parameters. For the FCC crystal with {nnn} =12,
{W}4 = -sqrt(14/143).(49/4096)/Pi^1.5 = -0.0006722136...
The optional keyword {wl/hat} will output the normalized third-order
invariants {Wlhat} (see Eq. 2.2 in "Steinhardt"_#Steinhardt)
for the same degrees as for the {Ql} parameters. For the FCC crystal
with {nnn} =12, {W}4hat = -7/3*sqrt(2/429) = -0.159317...The numerical
values of {Wlhat} for a range of commonly encountered high-symmetry
structures are given in Table I of "Steinhardt"_#Steinhardt, and these
can be reproduced with this keyword.
The optional keyword {components} will output the components of the The optional keyword {components} will output the components of the
normalized complex vector {Ybar_lm} of degree {ldegree}, which must be normalized complex vector {Ybar_lm} of degree {ldegree}, which must be
@ -82,7 +98,7 @@ particles, as discussed in "ten Wolde"_#tenWolde2.
The value of {Ql} is set to zero for atoms not in the The value of {Ql} is set to zero for atoms not in the
specified compute group, as well as for atoms that have less than specified compute group, as well as for atoms that have less than
{nnn} neighbors within the distance cutoff. {nnn} neighbors within the distance cutoff, unless {nnn} is NULL.
The neighbor list needed to compute this quantity is constructed each The neighbor list needed to compute this quantity is constructed each
time the calculation is performed (i.e. each time a snapshot of atoms time the calculation is performed (i.e. each time a snapshot of atoms
@ -108,6 +124,12 @@ This compute calculates a per-atom array with {nlvalues} columns,
giving the {Ql} values for each atom, which are real numbers on the giving the {Ql} values for each atom, which are real numbers on the
range 0 <= {Ql} <= 1. range 0 <= {Ql} <= 1.
If the keyword {wl} is set to yes, then the {Wl} values for each
atom will be added to the output array, which are real numbers.
If the keyword {wl/hat} is set to yes, then the {Wl_hat}
values for each atom will be added to the output array, which are real numbers.
If the keyword {components} is set, then the real and imaginary parts If the keyword {components} is set, then the real and imaginary parts
of each component of (normalized) {Ybar_lm} will be added to the of each component of (normalized) {Ybar_lm} will be added to the
output array in the following order: Re({Ybar_-m}) Im({Ybar_-m}) output array in the following order: Re({Ybar_-m}) Im({Ybar_-m})
@ -130,7 +152,8 @@ hexorder/atom"_compute_hexorder_atom.html
[Default:] [Default:]
The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, The option defaults are {cutoff} = pair style cutoff, {nnn} = 12,
{degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. {degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12,
{wl} = no, {wl/hat} = no, and {components} off
:line :line

91
examples/steinhardt/in.bcc Normal file
View File

@ -0,0 +1,91 @@
# Steinhardt-Nelson bond orientational order parameters for BCC
variable rcut equal 3.0
boundary p p p
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 5
# create geometry
lattice bcc 1.0
region box block 0 3 0 3 0 3
create_box 1 box
create_atoms 1 box
mass 1 1.0
# LJ potentials
pair_style lj/cut ${rcut}
pair_coeff * * 1.0 1.0 ${rcut}
# 14 neighbors, perfect crystal
compute qlwlhat all orientorder/atom degrees 6 2 4 6 8 10 12 nnn 14 wl/hat yes
compute avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] c_qlwlhat[6]
compute avwlhat all reduce ave c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] c_qlwlhat[11] c_qlwlhat[12]
thermo_style custom step temp epair etotal c_avql[*] c_avwlhat[*]
run 0
# check Q_l values
print " "
print "*******************************************************************"
print " "
print "Comparison with reference values of Q_l "
print " [Table I in W. Mickel, S. C. Kapfer,"
print " G. E. Schroeder-Turkand, K. Mecke, "
print " J. Chem. Phys. 138, 044501 (2013).]"
print " "
variable q2ref equal 0.0
variable q4ref equal 0.036
variable q6ref equal 0.511
variable q8ref equal 0.429
variable q10ref equal 0.195
variable q12ref equal 0.405
variable q2 equal c_avql[1]
variable q4 equal c_avql[2]
variable q6 equal c_avql[3]
variable q8 equal c_avql[4]
variable q10 equal c_avql[5]
variable q12 equal c_avql[6]
print "q2 = $(v_q2:%10.6f) delta = $(v_q2-v_q2ref:%10.4f)"
print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
# check W_l_hat values
print " "
print "Comparison with reference values of W_l_hat"
print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
print " Phys. Rev. B 28, 784 (1983).]"
print " "
variable w4hatref equal 0.159317
variable w6hatref equal 0.013161
variable w8hatref equal -0.058455
variable w10hatref equal -0.090130
variable w4hat equal c_avwlhat[2]
variable w6hat equal c_avwlhat[3]
variable w8hat equal c_avwlhat[4]
variable w10hat equal c_avwlhat[5]
print "w4hat = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)"
print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
print "w8hat = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)"
print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
print " "
print "*******************************************************************"
print " "

88
examples/steinhardt/in.fcc Normal file
View File

@ -0,0 +1,88 @@
# Steinhardt-Nelson bond orientational order parameters for FCC
variable rcut equal 3.0
boundary p p p
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 5
# create geometry
lattice fcc 1.0
region box block 0 3 0 3 0 3
create_box 1 box
create_atoms 1 box
mass 1 1.0
# LJ potentials
pair_style lj/cut ${rcut}
pair_coeff * * 1.0 1.0 ${rcut}
# 12 neighbors, perfect crystal
compute qlwlhat all orientorder/atom wl/hat yes
compute avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5]
compute avwlhat all reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10]
thermo_style custom step temp epair etotal c_avql[*] c_avwlhat[*]
run 0
# check Q_l values
print " "
print "*******************************************************************"
print " "
print "Comparison with reference values of Q_l "
print " [Table I in W. Mickel, S. C. Kapfer,"
print " G. E. Schroeder-Turkand, K. Mecke, "
print " J. Chem. Phys. 138, 044501 (2013).]"
print " "
variable q4ref equal 0.190
variable q6ref equal 0.575
variable q8ref equal 0.404
variable q10ref equal 0.013
variable q12ref equal 0.600
variable q4 equal c_avql[1]
variable q6 equal c_avql[2]
variable q8 equal c_avql[3]
variable q10 equal c_avql[4]
variable q12 equal c_avql[5]
print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
# check W_l_hat values
print " "
print "Comparison with reference values of W_l_hat"
print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
print " Phys. Rev. B 28, 784 (1983).]"
print " "
variable w4hatref equal -0.159316
variable w6hatref equal -0.013161
variable w8hatref equal 0.058454
variable w10hatref equal -0.090128
variable w4hat equal c_avwlhat[1]
variable w6hat equal c_avwlhat[2]
variable w8hat equal c_avwlhat[3]
variable w10hat equal c_avwlhat[4]
print "w4hat = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)"
print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
print "w8hat = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)"
print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
print " "
print "*******************************************************************"
print " "

106
examples/steinhardt/in.icos Normal file
View File

@ -0,0 +1,106 @@
# Steinhardt-Nelson bond orientational order parameters for icosahedral cluster
# W_6_hat is sensitive to icosohedral order
variable rcut equal 1.2 # a bit bigger than LJ Rmin
variable rcutred equal 0.75 # a bit bigger than 1/sqrt(2)
# create a perfect fcc crystallite
atom_style atomic
boundary s s s
lattice fcc 1.0 # neighbors at LJ Rmin
region box block 0 2 0 2 0 2
create_box 1 box
create_atoms 1 box
mass 1 1.0
region centralatom sphere 1 1 1 0.0 side in
group centralatom region centralatom
region mysphere sphere 1 1 1 ${rcutred} side out
delete_atoms region mysphere
# LJ potential
pair_style lj/cut 100.0
pair_coeff * * 1.0 1.0 100.0
# define output for central atom
compute qlwlhat all orientorder/atom wl/hat yes cutoff ${rcut} nnn NULL
compute avql centralatom reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5]
compute avwlhat centralatom reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10]
variable q6 equal c_avql[2]
variable w6hat equal c_avwlhat[2]
compute mype all pe/atom
compute centralatompe centralatom reduce ave c_mype
# gently equilibrate the crystallite
velocity all create 0.001 482748
fix 1 all nve
neighbor 0.3 bin
neigh_modify every 1 check no delay 0
timestep 0.003
thermo_style custom step temp epair etotal c_centralatompe v_q6 v_w6hat
thermo 10
run 10
# quench to icosehedral cluster
minimize 1.0e-10 1.0e-6 100 1000
# check Q_l values
print " "
print "*******************************************************************"
print " "
print "Comparison with reference values of Q_l "
print " [Table I in W. Mickel, S. C. Kapfer,"
print " G. E. Schroeder-Turkand, K. Mecke, "
print " J. Chem. Phys. 138, 044501 (2013).]"
print " "
variable q4ref equal 0.0
variable q6ref equal 0.663
variable q8ref equal 0.0
variable q10ref equal 0.363
variable q12ref equal 0.585
variable q4 equal c_avql[1]
variable q6 equal c_avql[2]
variable q8 equal c_avql[3]
variable q10 equal c_avql[4]
variable q12 equal c_avql[5]
print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
# check W_l_hat values
print " "
print "Comparison with reference values of W_l_hat"
print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
print " Phys. Rev. B 28, 784 (1983).]"
print " "
variable w6hatref equal -0.169754
variable w10hatref equal -0.093967
variable w4hat equal c_avwlhat[1]
variable w6hat equal c_avwlhat[2]
variable w8hat equal c_avwlhat[3]
variable w10hat equal c_avwlhat[4]
variable w12hat equal c_avwlhat[5]
print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
print " "
print "*******************************************************************"
print " "

459
src/compute_orientorder_atom.cpp
View File

@ -43,12 +43,14 @@ using namespace std;
#define MY_EPSILON (10.0*2.220446049250313e-16) #define MY_EPSILON (10.0*2.220446049250313e-16)
#endif #endif
#define QEPSILON 1.0e-6
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) : ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg),
qlist(NULL), distsq(NULL), nearest(NULL), rlist(NULL), qlist(NULL), distsq(NULL), nearest(NULL), rlist(NULL),
qnarray(NULL), qnm_r(NULL), qnm_i(NULL) qnarray(NULL), qnm_r(NULL), qnm_i(NULL), cglist(NULL)
{ {
if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command"); if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command");
@ -56,6 +58,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
nnn = 12; nnn = 12;
cutsq = 0.0; cutsq = 0.0;
wlflag = 0;
wlhatflag = 0;
qlcompflag = 0; qlcompflag = 0;
// specify which orders to request // specify which orders to request
@ -96,27 +100,39 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
if (iarg+nqlist > narg) if (iarg+nqlist > narg)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
qmax = 0; qmax = 0;
for (int iw = 0; iw < nqlist; iw++) { for (int il = 0; il < nqlist; il++) {
qlist[iw] = force->numeric(FLERR,arg[iarg+iw]); qlist[il] = force->numeric(FLERR,arg[iarg+il]);
if (qlist[iw] < 0) if (qlist[il] < 0)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
if (qlist[iw] > qmax) qmax = qlist[iw]; if (qlist[il] > qmax) qmax = qlist[il];
} }
iarg += nqlist; iarg += nqlist;
} else if (strcmp(arg[iarg],"wl") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command");
if (strcmp(arg[iarg+1],"yes") == 0) wlflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) wlflag = 0;
else error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"wl/hat") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command");
if (strcmp(arg[iarg+1],"yes") == 0) wlhatflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) wlhatflag = 0;
else error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"components") == 0) { } else if (strcmp(arg[iarg],"components") == 0) {
qlcompflag = 1; qlcompflag = 1;
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
qlcomp = force->numeric(FLERR,arg[iarg+1]); qlcomp = force->numeric(FLERR,arg[iarg+1]);
if (qlcomp <= 0)
error->all(FLERR,"Illegal compute orientorder/atom command");
iqlcomp = -1; iqlcomp = -1;
for (int iw = 0; iw < nqlist; iw++) for (int il = 0; il < nqlist; il++)
if (qlcomp == qlist[iw]) { if (qlcomp == qlist[il]) {
iqlcomp = iw; iqlcomp = il;
break; break;
} }
if (iqlcomp < 0) if (iqlcomp == -1)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"cutoff") == 0) { } else if (strcmp(arg[iarg],"cutoff") == 0) {
@ -130,8 +146,10 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
} else error->all(FLERR,"Illegal compute orientorder/atom command"); } else error->all(FLERR,"Illegal compute orientorder/atom command");
} }
if (qlcompflag) ncol = nqlist + 2*(2*qlcomp+1); ncol = nqlist;
else ncol = nqlist; if (wlflag) ncol += nqlist;
if (wlhatflag) ncol += nqlist;
if (qlcompflag) ncol += 2*(2*qlcomp+1);
peratom_flag = 1; peratom_flag = 1;
size_peratom_cols = ncol; size_peratom_cols = ncol;
@ -151,7 +169,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
memory->destroy(qlist); memory->destroy(qlist);
memory->destroy(qnm_r); memory->destroy(qnm_r);
memory->destroy(qnm_i); memory->destroy(qnm_i);
memory->destroy(cglist);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -166,8 +184,8 @@ void ComputeOrientOrderAtom::init()
error->all(FLERR,"Compute orientorder/atom cutoff is " error->all(FLERR,"Compute orientorder/atom cutoff is "
"longer than pairwise cutoff"); "longer than pairwise cutoff");
memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r"); memory->create(qnm_r,nqlist,2*qmax+1,"orientorder/atom:qnm_r");
memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i"); memory->create(qnm_i,nqlist,2*qmax+1,"orientorder/atom:qnm_i");
// need an occasional full neighbor list // need an occasional full neighbor list
@ -183,6 +201,8 @@ void ComputeOrientOrderAtom::init()
if (strcmp(modify->compute[i]->style,"orientorder/atom") == 0) count++; if (strcmp(modify->compute[i]->style,"orientorder/atom") == 0) count++;
if (count > 1 && comm->me == 0) if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute orientorder/atom"); error->warning(FLERR,"More than one compute orientorder/atom");
if (wlflag || wlhatflag) init_clebsch_gordan();
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -274,8 +294,8 @@ void ComputeOrientOrderAtom::compute_peratom()
// if not nnn neighbors, order parameter = 0; // if not nnn neighbors, order parameter = 0;
if ((ncount == 0) || (ncount < nnn)) { if ((ncount == 0) || (ncount < nnn)) {
for (int iw = 0; iw < nqlist; iw++) for (int jj = 0; jj < ncol; jj++)
qn[iw] = 0.0; qn[jj] = 0.0;
continue; continue;
} }
@ -287,6 +307,7 @@ void ComputeOrientOrderAtom::compute_peratom()
} }
calc_boop(rlist, ncount, qn, qlist, nqlist); calc_boop(rlist, ncount, qn, qlist, nqlist);
} }
} }
} }
@ -403,13 +424,12 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
void ComputeOrientOrderAtom::calc_boop(double **rlist, void ComputeOrientOrderAtom::calc_boop(double **rlist,
int ncount, double qn[], int ncount, double qn[],
int qlist[], int nqlist) { int qlist[], int nqlist) {
for (int iw = 0; iw < nqlist; iw++) {
int n = qlist[iw];
qn[iw] = 0.0; for (int il = 0; il < nqlist; il++) {
for(int m = 0; m < 2*n+1; m++) { int l = qlist[il];
qnm_r[iw][m] = 0.0; for(int m = 0; m < 2*l+1; m++) {
qnm_i[iw][m] = 0.0; qnm_r[il][m] = 0.0;
qnm_i[il][m] = 0.0;
} }
} }
@ -433,24 +453,24 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
expphi_i *= rxymaginv; expphi_i *= rxymaginv;
} }
for (int iw = 0; iw < nqlist; iw++) { for (int il = 0; il < nqlist; il++) {
int n = qlist[iw]; int l = qlist[il];
qnm_r[iw][n] += polar_prefactor(n, 0, costheta); qnm_r[il][l] += polar_prefactor(l, 0, costheta);
double expphim_r = expphi_r; double expphim_r = expphi_r;
double expphim_i = expphi_i; double expphim_i = expphi_i;
for(int m = 1; m <= +n; m++) { for(int m = 1; m <= +l; m++) {
double prefactor = polar_prefactor(n, m, costheta); double prefactor = polar_prefactor(l, m, costheta);
double c_r = prefactor * expphim_r; double c_r = prefactor * expphim_r;
double c_i = prefactor * expphim_i; double c_i = prefactor * expphim_i;
qnm_r[iw][m+n] += c_r; qnm_r[il][m+l] += c_r;
qnm_i[iw][m+n] += c_i; qnm_i[il][m+l] += c_i;
if(m & 1) { if(m & 1) {
qnm_r[iw][-m+n] -= c_r; qnm_r[il][-m+l] -= c_r;
qnm_i[iw][-m+n] += c_i; qnm_i[il][-m+l] += c_i;
} else { } else {
qnm_r[iw][-m+n] += c_r; qnm_r[il][-m+l] += c_r;
qnm_i[iw][-m+n] -= c_i; qnm_i[il][-m+l] -= c_i;
} }
double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i; double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r; double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
@ -461,30 +481,110 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
} }
} }
double fac = sqrt(MY_4PI) / ncount; // convert sums to averages
double normfac = 0.0;
for (int iw = 0; iw < nqlist; iw++) { double facn = 1.0 / ncount;
int n = qlist[iw]; for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for(int m = 0; m < 2*l+1; m++) {
qnm_r[il][m] *= facn;
qnm_i[il][m] *= facn;
}
}
// calculate Q_l
// 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++) {
int l = qlist[il];
double qnormfac = sqrt(MY_4PI/(2*l+1));
double qm_sum = 0.0; double qm_sum = 0.0;
for(int m = 0; m < 2*n+1; m++) { for(int m = 0; m < 2*l+1; m++)
qm_sum += qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]; qm_sum += qnm_r[il][m]*qnm_r[il][m] + qnm_i[il][m]*qnm_i[il][m];
// printf("Ylm^2 = %d %d %g\n",n,m, qn[jj++] = qnormfac * sqrt(qm_sum);
// qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]);
}
qn[iw] = fac * sqrt(qm_sum / (2*n+1));
if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum);
} }
// output of the complex vector // 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
// (Can try getting them from boop.py first)
// 5. [done]Compare to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop.py
// 6. [done]I get the right answer for W_l, but need to make sure that factor of 1/sqrt(l+1) is right for cglist
// 7. Add documentation
// 8. [done] run valgrind
// 9. [done] Add Wlhat
// 10. Update memory_usage()
// 11. Add exact FCC values for W_4, W_4_hat
// calculate W_l
if (wlflag) {
int idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
double wlsum = 0.0;
for(int m1 = 0; m1 < 2*l+1; m1++) {
for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
int m = m1 + m2 - l;
double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];
double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];
wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count];
idxcg_count++;
}
}
qn[jj++] = wlsum/sqrt(2*l+1);
}
}
// calculate W_l_hat
if (wlhatflag) {
int idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
double wlsum = 0.0;
for(int m1 = 0; m1 < 2*l+1; m1++) {
for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
int m = m1 + m2 - l;
double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];
double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];
wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count];
idxcg_count++;
}
}
// Whats = [w/(q/np.sqrt(np.pi * 4 / (2 * l + 1)))**3 if abs(q) > 1.0e-6 else 0.0 for l,q,w in zip(range(1,max_l+1),Qs,Ws)]
if (qn[il] < QEPSILON)
qn[jj++] = 0.0;
else {
double qnormfac = sqrt(MY_4PI/(2*l+1));
double qnfac = qnormfac/qn[il];
qn[jj++] = wlsum/sqrt(2*l+1)*(qnfac*qnfac*qnfac);
}
}
}
// Calculate components of Q_l, for l=qlcomp
if (qlcompflag) { if (qlcompflag) {
int j = nqlist; int il = iqlcomp;
for(int m = 0; m < 2*qlcomp+1; m++) { int l = qlcomp;
qn[j++] = qnm_r[iqlcomp][m] * normfac; if (qn[il] < QEPSILON)
qn[j++] = qnm_i[iqlcomp][m] * normfac; for(int m = 0; m < 2*l+1; m++) {
qn[jj++] = 0.0;
qn[jj++] = 0.0;
}
else {
double qnormfac = sqrt(MY_4PI/(2*l+1));
double qnfac = qnormfac/qn[il];
for(int m = 0; m < 2*l+1; m++) {
qn[jj++] = qnm_r[il][m] * qnfac;
qn[jj++] = qnm_i[il][m] * qnfac;
} }
} }
}
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -542,3 +642,258 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
return p; return p;
} }
/* ----------------------------------------------------------------------
assign Clebsch-Gordan coefficients
using the quasi-binomial formula VMK 8.2.1(3)
specialized for case j1=j2=j=l
------------------------------------------------------------------------- */
void ComputeOrientOrderAtom::init_clebsch_gordan()
{
double sum,dcg,sfaccg, sfac1, sfac2;
int m, aa2, bb2, cc2;
int ifac, idxcg_count;
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for(int m1 = 0; m1 < 2*l+1; m1++)
for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++)
idxcg_count++;
}
idxcg_max = idxcg_count;
memory->create(cglist, idxcg_max, "computeorientorderatom:cglist");
idxcg_count = 0;
for (int il = 0; il < nqlist; il++) {
int l = qlist[il];
for(int m1 = 0; m1 < 2*l+1; m1++) {
aa2 = m1 - l;
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++) {
ifac = z % 2 ? -1 : 1;
sum += ifac /
(factorial(z) *
factorial(l - z) *
factorial(l - aa2 - z) *
factorial(l + bb2 - z) *
factorial(aa2 + z) *
factorial(-bb2 + z));
}
cc2 = m - l;
sfaccg = sqrt(factorial(l + aa2) *
factorial(l - aa2) *
factorial(l + bb2) *
factorial(l - bb2) *
factorial(l + cc2) *
factorial(l - cc2) *
(2*l + 1));
sfac1 = factorial(3*l + 1);
sfac2 = factorial(l);
dcg = sqrt(sfac2*sfac2*sfac2 / sfac1);
cglist[idxcg_count] = sum * dcg * sfaccg;
idxcg_count++;
}
}
}
}
/* ----------------------------------------------------------------------
factorial n, wrapper for precomputed table
------------------------------------------------------------------------- */
double ComputeOrientOrderAtom::factorial(int n)
{
if (n < 0 || n > nmaxfactorial) {
char str[128];
sprintf(str, "Invalid argument to factorial %d", n);
error->all(FLERR, str);
}
return nfac_table[n];
}
/* ----------------------------------------------------------------------
factorial n table, size SNA::nmaxfactorial+1
------------------------------------------------------------------------- */
const double ComputeOrientOrderAtom::nfac_table[] = {
1,
1,
2,
6,
24,
120,
720,
5040,
40320,
362880,
3628800,
39916800,
479001600,
6227020800,
87178291200,
1307674368000,
20922789888000,
355687428096000,
6.402373705728e+15,
1.21645100408832e+17,
2.43290200817664e+18,
5.10909421717094e+19,
1.12400072777761e+21,
2.5852016738885e+22,
6.20448401733239e+23,
1.5511210043331e+25,
4.03291461126606e+26,
1.08888694504184e+28,
3.04888344611714e+29,
8.8417619937397e+30,
2.65252859812191e+32,
8.22283865417792e+33,
2.63130836933694e+35,
8.68331761881189e+36,
2.95232799039604e+38,
1.03331479663861e+40,
3.71993326789901e+41,
1.37637530912263e+43,
5.23022617466601e+44,
2.03978820811974e+46,
8.15915283247898e+47,
3.34525266131638e+49,
1.40500611775288e+51,
6.04152630633738e+52,
2.65827157478845e+54,
1.1962222086548e+56,
5.50262215981209e+57,
2.58623241511168e+59,
1.24139155925361e+61,
6.08281864034268e+62,
3.04140932017134e+64,
1.55111875328738e+66,
8.06581751709439e+67,
4.27488328406003e+69,
2.30843697339241e+71,
1.26964033536583e+73,
7.10998587804863e+74,
4.05269195048772e+76,
2.35056133128288e+78,
1.3868311854569e+80,
8.32098711274139e+81,
5.07580213877225e+83,
3.14699732603879e+85,
1.98260831540444e+87,
1.26886932185884e+89,
8.24765059208247e+90,
5.44344939077443e+92,
3.64711109181887e+94,
2.48003554243683e+96,
1.71122452428141e+98,
1.19785716699699e+100,
8.50478588567862e+101,
6.12344583768861e+103,
4.47011546151268e+105,
3.30788544151939e+107,
2.48091408113954e+109,
1.88549470166605e+111,
1.45183092028286e+113,
1.13242811782063e+115,
8.94618213078297e+116,
7.15694570462638e+118,
5.79712602074737e+120,
4.75364333701284e+122,
3.94552396972066e+124,
3.31424013456535e+126,
2.81710411438055e+128,
2.42270953836727e+130,
2.10775729837953e+132,
1.85482642257398e+134,
1.65079551609085e+136,
1.48571596448176e+138,
1.3520015276784e+140,
1.24384140546413e+142,
1.15677250708164e+144,
1.08736615665674e+146,
1.03299784882391e+148,
9.91677934870949e+149,
9.61927596824821e+151,
9.42689044888324e+153,
9.33262154439441e+155,
9.33262154439441e+157,
9.42594775983835e+159,
9.61446671503512e+161,
9.90290071648618e+163,
1.02990167451456e+166,
1.08139675824029e+168,
1.14628056373471e+170,
1.22652020319614e+172,
1.32464181945183e+174,
1.44385958320249e+176,
1.58824554152274e+178,
1.76295255109024e+180,
1.97450685722107e+182,
2.23119274865981e+184,
2.54355973347219e+186,
2.92509369349301e+188,
3.3931086844519e+190,
3.96993716080872e+192,
4.68452584975429e+194,
5.5745857612076e+196,
6.68950291344912e+198,
8.09429852527344e+200,
9.8750442008336e+202,
1.21463043670253e+205,
1.50614174151114e+207,
1.88267717688893e+209,
2.37217324288005e+211,
3.01266001845766e+213,
3.8562048236258e+215,
4.97450422247729e+217,
6.46685548922047e+219,
8.47158069087882e+221,
1.118248651196e+224,
1.48727070609069e+226,
1.99294274616152e+228,
2.69047270731805e+230,
3.65904288195255e+232,
5.01288874827499e+234,
6.91778647261949e+236,
9.61572319694109e+238,
1.34620124757175e+241,
1.89814375907617e+243,
2.69536413788816e+245,
3.85437071718007e+247,
5.5502938327393e+249,
8.04792605747199e+251,
1.17499720439091e+254,
1.72724589045464e+256,
2.55632391787286e+258,
3.80892263763057e+260,
5.71338395644585e+262,
8.62720977423323e+264,
1.31133588568345e+267,
2.00634390509568e+269,
3.08976961384735e+271,
4.78914290146339e+273,
7.47106292628289e+275,
1.17295687942641e+278,
1.85327186949373e+280,
2.94670227249504e+282,
4.71472363599206e+284,
7.59070505394721e+286,
1.22969421873945e+289,
2.0044015765453e+291,
3.28721858553429e+293,
5.42391066613159e+295,
9.00369170577843e+297,
1.503616514865e+300, // nmaxfactorial = 167
};

9
src/compute_orientorder_atom.h
View File

@ -33,7 +33,7 @@ class ComputeOrientOrderAtom : public Compute {
void compute_peratom(); void compute_peratom();
double memory_usage(); double memory_usage();
double cutsq; double cutsq;
int iqlcomp, qlcomp, qlcompflag; int iqlcomp, qlcomp, qlcompflag, wlflag, wlhatflag;
int *qlist; int *qlist;
int nqlist; int nqlist;
@ -55,6 +55,13 @@ class ComputeOrientOrderAtom : public Compute {
double polar_prefactor(int, int, double); double polar_prefactor(int, int, double);
double associated_legendre(int, int, double); double associated_legendre(int, int, double);
static const int nmaxfactorial = 167;
static const double nfac_table[];
double factorial(int);
void init_clebsch_gordan();
double *cglist; // Clebsch-Gordan coeffs
int idxcg_max;
}; };
} }