Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81064295
sna.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Sep 4, 01:11
Size
53 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 01:11 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20384830
Attached To
rLAMMPS lammps
sna.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Aidan Thompson, Christian Trott, SNL
------------------------------------------------------------------------- */
#include "sna.h"
#include <math.h>
#include "math_const.h"
#include "math_extra.h"
#include <string.h>
#include <stdlib.h>
#include "openmp_snap.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "atom.h"
using namespace std;
using namespace LAMMPS_NS;
using namespace MathConst;
/* ----------------------------------------------------------------------
this implementation is based on the method outlined
in Bartok[1], using formulae from VMK[2].
for the Clebsch-Gordan coefficients, we
convert the VMK half-integral labels
a, b, c, alpha, beta, gamma
to array offsets j1, j2, j, m1, m2, m
using the following relations:
j1 = 2*a
j2 = 2*b
j = 2*c
m1 = alpha+a 2*alpha = 2*m1 - j1
m2 = beta+b or 2*beta = 2*m2 - j2
m = gamma+c 2*gamma = 2*m - j
in this way:
-a <= alpha <= a
-b <= beta <= b
-c <= gamma <= c
becomes:
0 <= m1 <= j1
0 <= m2 <= j2
0 <= m <= j
and the requirement that
a+b+c be integral implies that
j1+j2+j must be even.
The requirement that:
gamma = alpha+beta
becomes:
2*m - j = 2*m1 - j1 + 2*m2 - j2
Similarly, for the Wigner U-functions U(J,m,m') we
convert the half-integral labels J,m,m' to
array offsets j,ma,mb:
j = 2*J
ma = J+m
mb = J+m'
so that:
0 <= j <= 2*Jmax
0 <= ma, mb <= j.
For the bispectrum components B(J1,J2,J) we convert to:
j1 = 2*J1
j2 = 2*J2
j = 2*J
and the requirement:
|J1-J2| <= J <= J1+J2, for j1+j2+j integral
becomes:
|j1-j2| <= j <= j1+j2, for j1+j2+j even integer
or
j = |j1-j2|, |j1-j2|+2,...,j1+j2-2,j1+j2
[1] Albert Bartok-Partay, "Gaussian Approximation..."
Doctoral Thesis, Cambrindge University, (2009)
[2] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,
"Quantum Theory of Angular Momentum," World Scientific (1988)
------------------------------------------------------------------------- */
SNA::SNA(LAMMPS* lmp, double rfac0_in,
int twojmax_in, int diagonalstyle_in, int use_shared_arrays_in,
double rmin0_in, int switch_flag_in) : Pointers(lmp)
{
wself = 1.0;
use_shared_arrays = use_shared_arrays_in;
rfac0 = rfac0_in;
rmin0 = rmin0_in;
switch_flag = switch_flag_in;
twojmax = twojmax_in;
diagonalstyle = diagonalstyle_in;
ncoeff = compute_ncoeff();
create_twojmax_arrays();
bvec = NULL;
dbvec = NULL;
memory->create(bvec, ncoeff, "pair:bvec");
memory->create(dbvec, ncoeff, 3, "pair:dbvec");
rij = NULL;
inside = NULL;
wj = NULL;
rcutij = NULL;
nmax = 0;
idxj = NULL;
#ifdef TIMING_INFO
timers = new double[20];
for(int i = 0; i < 20; i++) timers[i] = 0;
print = 0;
counter = 0;
#endif
build_indexlist();
}
/* ---------------------------------------------------------------------- */
SNA::~SNA()
{
if(!use_shared_arrays) {
destroy_twojmax_arrays();
memory->destroy(rij);
memory->destroy(inside);
memory->destroy(wj);
memory->destroy(rcutij);
memory->destroy(bvec);
memory->destroy(dbvec);
}
delete[] idxj;
}
void SNA::build_indexlist()
{
if(diagonalstyle == 0) {
int idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
idxj_count++;
// indexList can be changed here
idxj = new SNA_LOOPINDICES[idxj_count];
idxj_max = idxj_count;
idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
idxj[idxj_count].j1 = j1;
idxj[idxj_count].j2 = j2;
idxj[idxj_count].j = j;
idxj_count++;
}
}
if(diagonalstyle == 1) {
int idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j = 0; j <= MIN(twojmax, 2 * j1); j += 2) {
idxj_count++;
}
// indexList can be changed here
idxj = new SNA_LOOPINDICES[idxj_count];
idxj_max = idxj_count;
idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j = 0; j <= MIN(twojmax, 2 * j1); j += 2) {
idxj[idxj_count].j1 = j1;
idxj[idxj_count].j2 = j1;
idxj[idxj_count].j = j;
idxj_count++;
}
}
if(diagonalstyle == 2) {
int idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++) {
idxj_count++;
}
// indexList can be changed here
idxj = new SNA_LOOPINDICES[idxj_count];
idxj_max = idxj_count;
idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++) {
idxj[idxj_count].j1 = j1;
idxj[idxj_count].j2 = j1;
idxj[idxj_count].j = j1;
idxj_count++;
}
}
if(diagonalstyle == 3) {
int idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) idxj_count++;
// indexList can be changed here
idxj = new SNA_LOOPINDICES[idxj_count];
idxj_max = idxj_count;
idxj_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) {
idxj[idxj_count].j1 = j1;
idxj[idxj_count].j2 = j2;
idxj[idxj_count].j = j;
idxj_count++;
}
}
}
/* ---------------------------------------------------------------------- */
void SNA::init()
{
init_clebsch_gordan();
init_rootpqarray();
}
void SNA::grow_rij(int newnmax)
{
if(newnmax <= nmax) return;
nmax = newnmax;
if(!use_shared_arrays) {
memory->destroy(rij);
memory->destroy(inside);
memory->destroy(wj);
memory->destroy(rcutij);
memory->create(rij, nmax, 3, "pair:rij");
memory->create(inside, nmax, "pair:inside");
memory->create(wj, nmax, "pair:wj");
memory->create(rcutij, nmax, "pair:rcutij");
}
}
/* ----------------------------------------------------------------------
compute Ui by summing over neighbors j
------------------------------------------------------------------------- */
void SNA::compute_ui(int jnum)
{
double rsq, r, x, y, z, z0, theta0;
// utot(j,ma,mb) = 0 for all j,ma,ma
// utot(j,ma,ma) = 1 for all j,ma
// for j in neighbors of i:
// compute r0 = (x,y,z,z0)
// utot(j,ma,mb) += u(r0;j,ma,mb) for all j,ma,mb
zero_uarraytot();
addself_uarraytot(wself);
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
for(int j = 0; j < jnum; j++) {
x = rij[j][0];
y = rij[j][1];
z = rij[j][2];
rsq = x * x + y * y + z * z;
r = sqrt(rsq);
theta0 = (r - rmin0) * rfac0 * MY_PI / (rcutij[j] - rmin0);
// theta0 = (r - rmin0) * rscale0;
z0 = r / tan(theta0);
compute_uarray(x, y, z, z0, r);
add_uarraytot(r, wj[j], rcutij[j]);
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[0] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
void SNA::compute_ui_omp(int jnum, int sub_threads)
{
double rsq, r, x, y, z, z0, theta0;
// utot(j,ma,mb) = 0 for all j,ma,ma
// utot(j,ma,ma) = 1 for all j,ma
// for j in neighbors of i:
// compute r0 = (x,y,z,z0)
// utot(j,ma,mb) += u(r0;j,ma,mb) for all j,ma,mb
zero_uarraytot();
addself_uarraytot(wself);
for(int j = 0; j < jnum; j++) {
x = rij[j][0];
y = rij[j][1];
z = rij[j][2];
rsq = x * x + y * y + z * z;
r = sqrt(rsq);
theta0 = (r - rmin0) * rfac0 * MY_PI / (rcutij[j] - rmin0);
// theta0 = (r - rmin0) * rscale0;
z0 = r / tan(theta0);
omp_set_num_threads(sub_threads);
#if defined(_OPENMP)
#pragma omp parallel shared(x,y,z,z0,r,sub_threads) default(none)
#endif
{
compute_uarray_omp(x, y, z, z0, r, sub_threads);
}
add_uarraytot(r, wj[j], rcutij[j]);
}
}
/* ----------------------------------------------------------------------
compute Zi by summing over products of Ui
------------------------------------------------------------------------- */
void SNA::compute_zi()
{
// for j1 = 0,...,twojmax
// for j2 = 0,twojmax
// for j = |j1-j2|,Min(twojmax,j1+j2),2
// for ma = 0,...,j
// for mb = 0,...,jmid
// z(j1,j2,j,ma,mb) = 0
// for ma1 = Max(0,ma+(j1-j2-j)/2),Min(j1,ma+(j1+j2-j)/2)
// sumb1 = 0
// ma2 = ma-ma1+(j1+j2-j)/2;
// for mb1 = Max(0,mb+(j1-j2-j)/2),Min(j1,mb+(j1+j2-j)/2)
// mb2 = mb-mb1+(j1+j2-j)/2;
// sumb1 += cg(j1,mb1,j2,mb2,j) *
// u(j1,ma1,mb1) * u(j2,ma2,mb2)
// z(j1,j2,j,ma,mb) += sumb1*cg(j1,ma1,j2,ma2,j)
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
// compute_dbidrj() requires full j1/j2/j chunk of z elements
// use zarray j1/j2 symmetry
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) {
double sumb1_r, sumb1_i;
int ma2, mb2;
for(int mb = 0; 2*mb <= j; mb++)
for(int ma = 0; ma <= j; ma++) {
zarray_r[j1][j2][j][ma][mb] = 0.0;
zarray_i[j1][j2][j][ma][mb] = 0.0;
for(int ma1 = MAX(0, (2 * ma - j - j2 + j1) / 2);
ma1 <= MIN(j1, (2 * ma - j + j2 + j1) / 2); ma1++) {
sumb1_r = 0.0;
sumb1_i = 0.0;
ma2 = (2 * ma - j - (2 * ma1 - j1) + j2) / 2;
for(int mb1 = MAX(0, (2 * mb - j - j2 + j1) / 2);
mb1 <= MIN(j1, (2 * mb - j + j2 + j1) / 2); mb1++) {
mb2 = (2 * mb - j - (2 * mb1 - j1) + j2) / 2;
sumb1_r += cgarray[j1][j2][j][mb1][mb2] *
(uarraytot_r[j1][ma1][mb1] * uarraytot_r[j2][ma2][mb2] -
uarraytot_i[j1][ma1][mb1] * uarraytot_i[j2][ma2][mb2]);
sumb1_i += cgarray[j1][j2][j][mb1][mb2] *
(uarraytot_r[j1][ma1][mb1] * uarraytot_i[j2][ma2][mb2] +
uarraytot_i[j1][ma1][mb1] * uarraytot_r[j2][ma2][mb2]);
} // end loop over mb1
zarray_r[j1][j2][j][ma][mb] +=
sumb1_r * cgarray[j1][j2][j][ma1][ma2];
zarray_i[j1][j2][j][ma][mb] +=
sumb1_i * cgarray[j1][j2][j][ma1][ma2];
} // end loop over ma1
} // end loop over ma, mb
} // end loop over j
} // end loop over j1, j2
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[1] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
void SNA::compute_zi_omp(int sub_threads)
{
// for j1 = 0,...,twojmax
// for j2 = 0,twojmax
// for j = |j1-j2|,Min(twojmax,j1+j2),2
// for ma = 0,...,j
// for mb = 0,...,j
// z(j1,j2,j,ma,mb) = 0
// for ma1 = Max(0,ma+(j1-j2-j)/2),Min(j1,ma+(j1+j2-j)/2)
// sumb1 = 0
// ma2 = ma-ma1+(j1+j2-j)/2;
// for mb1 = Max(0,mb+(j1-j2-j)/2),Min(j1,mb+(j1+j2-j)/2)
// mb2 = mb-mb1+(j1+j2-j)/2;
// sumb1 += cg(j1,mb1,j2,mb2,j) *
// u(j1,ma1,mb1) * u(j2,ma2,mb2)
// z(j1,j2,j,ma,mb) += sumb1*cg(j1,ma1,j2,ma2,j)
if(omp_in_parallel())
omp_set_num_threads(sub_threads);
// compute_dbidrj() requires full j1/j2/j chunk of z elements
// use zarray j1/j2 symmetry
#if defined(_OPENMP)
#pragma omp parallel for schedule(auto) default(none)
#endif
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
double sumb1_r, sumb1_i;
int ma2, mb2;
for(int ma = 0; ma <= j; ma++)
for(int mb = 0; mb <= j; mb++) {
zarray_r[j1][j2][j][ma][mb] = 0.0;
zarray_i[j1][j2][j][ma][mb] = 0.0;
for(int ma1 = MAX(0, (2 * ma - j - j2 + j1) / 2);
ma1 <= MIN(j1, (2 * ma - j + j2 + j1) / 2); ma1++) {
sumb1_r = 0.0;
sumb1_i = 0.0;
ma2 = (2 * ma - j - (2 * ma1 - j1) + j2) / 2;
for(int mb1 = MAX(0, (2 * mb - j - j2 + j1) / 2);
mb1 <= MIN(j1, (2 * mb - j + j2 + j1) / 2); mb1++) {
mb2 = (2 * mb - j - (2 * mb1 - j1) + j2) / 2;
sumb1_r += cgarray[j1][j2][j][mb1][mb2] *
(uarraytot_r[j1][ma1][mb1] * uarraytot_r[j2][ma2][mb2] -
uarraytot_i[j1][ma1][mb1] * uarraytot_i[j2][ma2][mb2]);
sumb1_i += cgarray[j1][j2][j][mb1][mb2] *
(uarraytot_r[j1][ma1][mb1] * uarraytot_i[j2][ma2][mb2] +
uarraytot_i[j1][ma1][mb1] * uarraytot_r[j2][ma2][mb2]);
}
zarray_r[j1][j2][j][ma][mb] +=
sumb1_r * cgarray[j1][j2][j][ma1][ma2];
zarray_i[j1][j2][j][ma][mb] +=
sumb1_i * cgarray[j1][j2][j][ma1][ma2];
}
}
}
}
/* ----------------------------------------------------------------------
compute Bi by summing conj(Ui)*Zi
------------------------------------------------------------------------- */
void SNA::compute_bi()
{
// for j1 = 0,...,twojmax
// for j2 = 0,twojmax
// for j = |j1-j2|,Min(twojmax,j1+j2),2
// b(j1,j2,j) = 0
// for mb = 0,...,jmid
// for ma = 0,...,j
// b(j1,j2,j) +=
// 2*Conj(u(j,ma,mb))*z(j1,j2,j,ma,mb)
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++) {
for(int j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2) {
barray[j1][j2][j] = 0.0;
for(int mb = 0; 2*mb < j; mb++) {
for(int ma = 0; ma <= j; ma++) {
barray[j1][j2][j] +=
uarraytot_r[j][ma][mb] * zarray_r[j1][j2][j][ma][mb] +
uarraytot_i[j][ma][mb] * zarray_i[j1][j2][j][ma][mb];
}
}
// For j even, special treatment for middle column
if (j%2 == 0) {
int mb = j/2;
for(int ma = 0; ma < mb; ma++)
barray[j1][j2][j] +=
uarraytot_r[j][ma][mb] * zarray_r[j1][j2][j][ma][mb] +
uarraytot_i[j][ma][mb] * zarray_i[j1][j2][j][ma][mb];
int ma = mb;
barray[j1][j2][j] +=
(uarraytot_r[j][ma][mb] * zarray_r[j1][j2][j][ma][mb] +
uarraytot_i[j][ma][mb] * zarray_i[j1][j2][j][ma][mb])*0.5;
}
barray[j1][j2][j] *= 2.0;
}
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[2] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
/* ----------------------------------------------------------------------
copy Bi array to a vector
------------------------------------------------------------------------- */
void SNA::copy_bi2bvec()
{
int ncount, j1, j2, j;
ncount = 0;
for(j1 = 0; j1 <= twojmax; j1++)
if(diagonalstyle == 0) {
for(j2 = 0; j2 <= j1; j2++)
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2) {
bvec[ncount] = barray[j1][j2][j];
ncount++;
}
} else if(diagonalstyle == 1) {
j2 = j1;
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2) {
bvec[ncount] = barray[j1][j2][j];
ncount++;
}
} else if(diagonalstyle == 2) {
j = j2 = j1;
bvec[ncount] = barray[j1][j2][j];
ncount++;
} else if(diagonalstyle == 3) {
for(j2 = 0; j2 <= j1; j2++)
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) {
bvec[ncount] = barray[j1][j2][j];
ncount++;
}
}
}
/* ----------------------------------------------------------------------
calculate derivative of Ui w.r.t. atom j
------------------------------------------------------------------------- */
void SNA::compute_duidrj(double* rij, double wj, double rcut)
{
double rsq, r, x, y, z, z0, theta0, cs, sn;
double dz0dr;
x = rij[0];
y = rij[1];
z = rij[2];
rsq = x * x + y * y + z * z;
r = sqrt(rsq);
double rscale0 = rfac0 * MY_PI / (rcut - rmin0);
theta0 = (r - rmin0) * rscale0;
cs = cos(theta0);
sn = sin(theta0);
z0 = r * cs / sn;
dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq;
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut);
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[3] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
/* ----------------------------------------------------------------------
calculate derivative of Bi w.r.t. atom j
variant using indexlist for j1,j2,j
variant not using symmetry relation
------------------------------------------------------------------------- */
void SNA::compute_dbidrj_nonsymm()
{
// for j1 = 0,...,twojmax
// for j2 = 0,twojmax
// for j = |j1-j2|,Min(twojmax,j1+j2),2
// dbdr(j1,j2,j) = 0
// for ma = 0,...,j
// for mb = 0,...,j
// dzdr = 0
// for ma1 = Max(0,ma+(j1-j2-j)/2),Min(j1,ma+(j1+j2-j)/2)
// sumb1 = 0
// ma2 = ma-ma1+(j1+j2-j)/2;
// for mb1 = Max(0,mb+(j1-j2-j)/2),Min(j1,mb+(j1+j2-j)/2)
// mb2 = mb-mb1+(j1+j2-j)/2;
// sumb1 += cg(j1,mb1,j2,mb2,j) *
// (dudr(j1,ma1,mb1) * u(j2,ma2,mb2) +
// u(j1,ma1,mb1) * dudr(j2,ma2,mb2))
// dzdr += sumb1*cg(j1,ma1,j2,ma2,j)
// dbdr(j1,j2,j) +=
// Conj(dudr(j,ma,mb))*z(j1,j2,j,ma,mb) +
// Conj(u(j,ma,mb))*dzdr
double* dbdr;
double* dudr_r, *dudr_i;
double sumb1_r[3], sumb1_i[3], dzdr_r[3], dzdr_i[3];
int ma2;
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
for(int JJ = 0; JJ < idxj_max; JJ++) {
const int j1 = idxj[JJ].j1;
const int j2 = idxj[JJ].j2;
const int j = idxj[JJ].j;
dbdr = dbarray[j1][j2][j];
dbdr[0] = 0.0;
dbdr[1] = 0.0;
dbdr[2] = 0.0;
double** *j1duarray_r = duarray_r[j1];
double** *j2duarray_r = duarray_r[j2];
double** *j1duarray_i = duarray_i[j1];
double** *j2duarray_i = duarray_i[j2];
double** j1uarraytot_r = uarraytot_r[j1];
double** j2uarraytot_r = uarraytot_r[j2];
double** j1uarraytot_i = uarraytot_i[j1];
double** j2uarraytot_i = uarraytot_i[j2];
double** j1j2jcgarray = cgarray[j1][j2][j];
for(int ma = 0; ma <= j; ma++)
for(int mb = 0; mb <= j; mb++) {
dzdr_r[0] = 0.0;
dzdr_r[1] = 0.0;
dzdr_r[2] = 0.0;
dzdr_i[0] = 0.0;
dzdr_i[1] = 0.0;
dzdr_i[2] = 0.0;
const int max_mb1 = MIN(j1, (2 * mb - j + j2 + j1) / 2) + 1;
const int max_ma1 = MIN(j1, (2 * ma - j + j2 + j1) / 2) + 1;
for(int ma1 = MAX(0, (2 * ma - j - j2 + j1) / 2);
ma1 < max_ma1; ma1++) {
ma2 = (2 * ma - j - (2 * ma1 - j1) + j2) / 2;
sumb1_r[0] = 0.0;
sumb1_r[1] = 0.0;
sumb1_r[2] = 0.0;
sumb1_i[0] = 0.0;
sumb1_i[1] = 0.0;
sumb1_i[2] = 0.0;
//inside loop 54 operations (mul and add)
for(int mb1 = MAX(0, (2 * mb - j - j2 + j1) / 2),
mb2 = mb + (j1 + j2 - j) / 2 - mb1;
mb1 < max_mb1; mb1++, mb2--) {
double* dudr1_r, *dudr1_i, *dudr2_r, *dudr2_i;
dudr1_r = j1duarray_r[ma1][mb1];
dudr2_r = j2duarray_r[ma2][mb2];
dudr1_i = j1duarray_i[ma1][mb1];
dudr2_i = j2duarray_i[ma2][mb2];
const double cga_mb1mb2 = j1j2jcgarray[mb1][mb2];
const double uat_r_ma2mb2 = cga_mb1mb2 * j2uarraytot_r[ma2][mb2];
const double uat_r_ma1mb1 = cga_mb1mb2 * j1uarraytot_r[ma1][mb1];
const double uat_i_ma2mb2 = cga_mb1mb2 * j2uarraytot_i[ma2][mb2];
const double uat_i_ma1mb1 = cga_mb1mb2 * j1uarraytot_i[ma1][mb1];
for(int k = 0; k < 3; k++) {
sumb1_r[k] += dudr1_r[k] * uat_r_ma2mb2;
sumb1_r[k] -= dudr1_i[k] * uat_i_ma2mb2;
sumb1_i[k] += dudr1_r[k] * uat_i_ma2mb2;
sumb1_i[k] += dudr1_i[k] * uat_r_ma2mb2;
sumb1_r[k] += dudr2_r[k] * uat_r_ma1mb1;
sumb1_r[k] -= dudr2_i[k] * uat_i_ma1mb1;
sumb1_i[k] += dudr2_r[k] * uat_i_ma1mb1;
sumb1_i[k] += dudr2_i[k] * uat_r_ma1mb1;
}
} // end loop over mb1,mb2
// dzdr += sumb1*cg(j1,ma1,j2,ma2,j)
dzdr_r[0] += sumb1_r[0] * j1j2jcgarray[ma1][ma2];
dzdr_r[1] += sumb1_r[1] * j1j2jcgarray[ma1][ma2];
dzdr_r[2] += sumb1_r[2] * j1j2jcgarray[ma1][ma2];
dzdr_i[0] += sumb1_i[0] * j1j2jcgarray[ma1][ma2];
dzdr_i[1] += sumb1_i[1] * j1j2jcgarray[ma1][ma2];
dzdr_i[2] += sumb1_i[2] * j1j2jcgarray[ma1][ma2];
} // end loop over ma1,ma2
// dbdr(j1,j2,j) +=
// Conj(dudr(j,ma,mb))*z(j1,j2,j,ma,mb) +
// Conj(u(j,ma,mb))*dzdr
dudr_r = duarray_r[j][ma][mb];
dudr_i = duarray_i[j][ma][mb];
for(int k = 0; k < 3; k++)
dbdr[k] +=
(dudr_r[k] * zarray_r[j1][j2][j][ma][mb] +
dudr_i[k] * zarray_i[j1][j2][j][ma][mb]) +
(uarraytot_r[j][ma][mb] * dzdr_r[k] +
uarraytot_i[j][ma][mb] * dzdr_i[k]);
} //end loop over ma mb
} //end loop over j1 j2 j
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[4] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
/* ----------------------------------------------------------------------
calculate derivative of Bi w.r.t. atom j
variant using indexlist for j1,j2,j
variant using symmetry relation
------------------------------------------------------------------------- */
void SNA::compute_dbidrj()
{
// for j1 = 0,...,twojmax
// for j2 = 0,twojmax
// for j = |j1-j2|,Min(twojmax,j1+j2),2
// zdb = 0
// for mb = 0,...,jmid
// for ma = 0,...,j
// zdb +=
// Conj(dudr(j,ma,mb))*z(j1,j2,j,ma,mb)
// dbdr(j1,j2,j) += 2*zdb
// zdb = 0
// for mb1 = 0,...,j1mid
// for ma1 = 0,...,j1
// zdb +=
// Conj(dudr(j1,ma1,mb1))*z(j1,j2,j,ma1,mb1)
// dbdr(j1,j2,j) += 2*zdb*(j+1)/(j1+1)
// zdb = 0
// for mb2 = 0,...,j2mid
// for ma2 = 0,...,j2
// zdb +=
// Conj(dudr(j2,ma2,mb2))*z(j1,j,j2,ma2,mb2)
// dbdr(j1,j2,j) += 2*zdb*(j+1)/(j2+1)
double* dbdr;
double* dudr_r, *dudr_i;
double sumzdu_r[3];
double** jjjzarray_r;
double** jjjzarray_i;
double jjjmambzarray_r;
double jjjmambzarray_i;
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
#endif
for(int JJ = 0; JJ < idxj_max; JJ++) {
const int j1 = idxj[JJ].j1;
const int j2 = idxj[JJ].j2;
const int j = idxj[JJ].j;
dbdr = dbarray[j1][j2][j];
dbdr[0] = 0.0;
dbdr[1] = 0.0;
dbdr[2] = 0.0;
// Sum terms Conj(dudr(j,ma,mb))*z(j1,j2,j,ma,mb)
for(int k = 0; k < 3; k++)
sumzdu_r[k] = 0.0;
// use zarray j1/j2 symmetry (optional)
if (j1 >= j2) {
jjjzarray_r = zarray_r[j1][j2][j];
jjjzarray_i = zarray_i[j1][j2][j];
} else {
jjjzarray_r = zarray_r[j2][j1][j];
jjjzarray_i = zarray_i[j2][j1][j];
}
for(int mb = 0; 2*mb < j; mb++)
for(int ma = 0; ma <= j; ma++) {
dudr_r = duarray_r[j][ma][mb];
dudr_i = duarray_i[j][ma][mb];
jjjmambzarray_r = jjjzarray_r[ma][mb];
jjjmambzarray_i = jjjzarray_i[ma][mb];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
} //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++) {
dudr_r = duarray_r[j][ma][mb];
dudr_i = duarray_i[j][ma][mb];
jjjmambzarray_r = jjjzarray_r[ma][mb];
jjjmambzarray_i = jjjzarray_i[ma][mb];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
}
int ma = mb;
dudr_r = duarray_r[j][ma][mb];
dudr_i = duarray_i[j][ma][mb];
jjjmambzarray_r = jjjzarray_r[ma][mb];
jjjmambzarray_i = jjjzarray_i[ma][mb];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
(dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i)*0.5;
} // end if jeven
for(int k = 0; k < 3; k++)
dbdr[k] += 2.0*sumzdu_r[k];
// Sum over Conj(dudr(j1,ma1,mb1))*z(j,j2,j1,ma1,mb1)
double j1fac = (j+1)/(j1+1.0);
for(int k = 0; k < 3; k++)
sumzdu_r[k] = 0.0;
// use zarray j1/j2 symmetry (optional)
if (j >= j2) {
jjjzarray_r = zarray_r[j][j2][j1];
jjjzarray_i = zarray_i[j][j2][j1];
} else {
jjjzarray_r = zarray_r[j2][j][j1];
jjjzarray_i = zarray_i[j2][j][j1];
}
for(int mb1 = 0; 2*mb1 < j1; mb1++)
for(int ma1 = 0; ma1 <= j1; ma1++) {
dudr_r = duarray_r[j1][ma1][mb1];
dudr_i = duarray_i[j1][ma1][mb1];
jjjmambzarray_r = jjjzarray_r[ma1][mb1];
jjjmambzarray_i = jjjzarray_i[ma1][mb1];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
} //end loop over ma1 mb1
// For j1 even, handle middle column
if (j1%2 == 0) {
int mb1 = j1/2;
for(int ma1 = 0; ma1 < mb1; ma1++) {
dudr_r = duarray_r[j1][ma1][mb1];
dudr_i = duarray_i[j1][ma1][mb1];
jjjmambzarray_r = jjjzarray_r[ma1][mb1];
jjjmambzarray_i = jjjzarray_i[ma1][mb1];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
}
int ma1 = mb1;
dudr_r = duarray_r[j1][ma1][mb1];
dudr_i = duarray_i[j1][ma1][mb1];
jjjmambzarray_r = jjjzarray_r[ma1][mb1];
jjjmambzarray_i = jjjzarray_i[ma1][mb1];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
(dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i)*0.5;
} // end if j1even
for(int k = 0; k < 3; k++)
dbdr[k] += 2.0*sumzdu_r[k]*j1fac;
// Sum over Conj(dudr(j2,ma2,mb2))*z(j1,j,j2,ma2,mb2)
double j2fac = (j+1)/(j2+1.0);
for(int k = 0; k < 3; k++)
sumzdu_r[k] = 0.0;
// use zarray j1/j2 symmetry (optional)
if (j1 >= j) {
jjjzarray_r = zarray_r[j1][j][j2];
jjjzarray_i = zarray_i[j1][j][j2];
} else {
jjjzarray_r = zarray_r[j][j1][j2];
jjjzarray_i = zarray_i[j][j1][j2];
}
for(int mb2 = 0; 2*mb2 < j2; mb2++)
for(int ma2 = 0; ma2 <= j2; ma2++) {
dudr_r = duarray_r[j2][ma2][mb2];
dudr_i = duarray_i[j2][ma2][mb2];
jjjmambzarray_r = jjjzarray_r[ma2][mb2];
jjjmambzarray_i = jjjzarray_i[ma2][mb2];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
} //end loop over ma2 mb2
// For j2 even, handle middle column
if (j2%2 == 0) {
int mb2 = j2/2;
for(int ma2 = 0; ma2 < mb2; ma2++) {
dudr_r = duarray_r[j2][ma2][mb2];
dudr_i = duarray_i[j2][ma2][mb2];
jjjmambzarray_r = jjjzarray_r[ma2][mb2];
jjjmambzarray_i = jjjzarray_i[ma2][mb2];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i;
}
int ma2 = mb2;
dudr_r = duarray_r[j2][ma2][mb2];
dudr_i = duarray_i[j2][ma2][mb2];
jjjmambzarray_r = jjjzarray_r[ma2][mb2];
jjjmambzarray_i = jjjzarray_i[ma2][mb2];
for(int k = 0; k < 3; k++)
sumzdu_r[k] +=
(dudr_r[k] * jjjmambzarray_r +
dudr_i[k] * jjjmambzarray_i)*0.5;
} // end if j2even
for(int k = 0; k < 3; k++)
dbdr[k] += 2.0*sumzdu_r[k]*j2fac;
} //end loop over j1 j2 j
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &endtime);
timers[4] += (endtime.tv_sec - starttime.tv_sec + 1.0 *
(endtime.tv_nsec - starttime.tv_nsec) / 1000000000);
#endif
}
/* ----------------------------------------------------------------------
copy Bi derivatives into a vector
------------------------------------------------------------------------- */
void SNA::copy_dbi2dbvec()
{
int ncount, j1, j2, j;
ncount = 0;
for(j1 = 0; j1 <= twojmax; j1++) {
if(diagonalstyle == 0) {
for(j2 = 0; j2 <= j1; j2++)
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2) {
dbvec[ncount][0] = dbarray[j1][j2][j][0];
dbvec[ncount][1] = dbarray[j1][j2][j][1];
dbvec[ncount][2] = dbarray[j1][j2][j][2];
ncount++;
}
} else if(diagonalstyle == 1) {
j2 = j1;
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2) {
dbvec[ncount][0] = dbarray[j1][j2][j][0];
dbvec[ncount][1] = dbarray[j1][j2][j][1];
dbvec[ncount][2] = dbarray[j1][j2][j][2];
ncount++;
}
} else if(diagonalstyle == 2) {
j = j2 = j1;
dbvec[ncount][0] = dbarray[j1][j2][j][0];
dbvec[ncount][1] = dbarray[j1][j2][j][1];
dbvec[ncount][2] = dbarray[j1][j2][j][2];
ncount++;
} else if(diagonalstyle == 3) {
for(j2 = 0; j2 <= j1; j2++)
for(j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) {
dbvec[ncount][0] = dbarray[j1][j2][j][0];
dbvec[ncount][1] = dbarray[j1][j2][j][1];
dbvec[ncount][2] = dbarray[j1][j2][j][2];
ncount++;
}
}
}
}
/* ---------------------------------------------------------------------- */
void SNA::zero_uarraytot()
{
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++)
for (int mb = 0; mb <= j; mb++) {
uarraytot_r[j][ma][mb] = 0.0;
uarraytot_i[j][ma][mb] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
void SNA::addself_uarraytot(double wself_in)
{
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++) {
uarraytot_r[j][ma][ma] = wself_in;
uarraytot_i[j][ma][ma] = 0.0;
}
}
/* ----------------------------------------------------------------------
add Wigner U-functions for one neighbor to the total
------------------------------------------------------------------------- */
void SNA::add_uarraytot(double r, double wj, double rcut)
{
double sfac;
sfac = compute_sfac(r, rcut);
sfac *= wj;
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++)
for (int mb = 0; mb <= j; mb++) {
uarraytot_r[j][ma][mb] +=
sfac * uarray_r[j][ma][mb];
uarraytot_i[j][ma][mb] +=
sfac * uarray_i[j][ma][mb];
}
}
void SNA::add_uarraytot_omp(double r, double wj, double rcut)
{
double sfac;
sfac = compute_sfac(r, rcut);
sfac *= wj;
#if defined(_OPENMP)
#pragma omp for
#endif
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++)
for (int mb = 0; mb <= j; mb++) {
uarraytot_r[j][ma][mb] +=
sfac * uarray_r[j][ma][mb];
uarraytot_i[j][ma][mb] +=
sfac * uarray_i[j][ma][mb];
}
}
/* ----------------------------------------------------------------------
compute Wigner U-functions for one neighbor
------------------------------------------------------------------------- */
void SNA::compute_uarray(double x, double y, double z,
double z0, double r)
{
double r0inv;
double a_r, b_r, a_i, b_i;
double rootpq;
// compute Cayley-Klein parameters for unit quaternion
r0inv = 1.0 / sqrt(r * r + z0 * z0);
a_r = r0inv * z0;
a_i = -r0inv * z;
b_r = r0inv * y;
b_i = -r0inv * x;
// VMK Section 4.8.2
uarray_r[0][0][0] = 1.0;
uarray_i[0][0][0] = 0.0;
for (int j = 1; j <= twojmax; j++) {
// fill in left side of matrix layer from previous layer
for (int mb = 0; 2*mb <= j; mb++) {
uarray_r[j][0][mb] = 0.0;
uarray_i[j][0][mb] = 0.0;
for (int ma = 0; ma < j; ma++) {
rootpq = rootpqarray[j - ma][j - mb];
uarray_r[j][ma][mb] +=
rootpq *
(a_r * uarray_r[j - 1][ma][mb] +
a_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma][mb] +=
rootpq *
(a_r * uarray_i[j - 1][ma][mb] -
a_i * uarray_r[j - 1][ma][mb]);
rootpq = rootpqarray[ma + 1][j - mb];
uarray_r[j][ma + 1][mb] =
-rootpq *
(b_r * uarray_r[j - 1][ma][mb] +
b_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma + 1][mb] =
-rootpq *
(b_r * uarray_i[j - 1][ma][mb] -
b_i * uarray_r[j - 1][ma][mb]);
}
}
// copy left side to right side with inversion symmetry VMK 4.4(2)
// u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])
int mbpar = -1;
for (int mb = 0; 2*mb <= j; mb++) {
mbpar = -mbpar;
int mapar = -mbpar;
for (int ma = 0; ma <= j; ma++) {
mapar = -mapar;
if (mapar == 1) {
uarray_r[j][j-ma][j-mb] = uarray_r[j][ma][mb];
uarray_i[j][j-ma][j-mb] = -uarray_i[j][ma][mb];
} else {
uarray_r[j][j-ma][j-mb] = -uarray_r[j][ma][mb];
uarray_i[j][j-ma][j-mb] = uarray_i[j][ma][mb];
}
}
}
}
}
void SNA::compute_uarray_omp(double x, double y, double z,
double z0, double r, int sub_threads)
{
double r0inv;
double a_r, b_r, a_i, b_i;
double rootpq;
// compute Cayley-Klein parameters for unit quaternion
r0inv = 1.0 / sqrt(r * r + z0 * z0);
a_r = r0inv * z0;
a_i = -r0inv * z;
b_r = r0inv * y;
b_i = -r0inv * x;
// VMK Section 4.8.2
uarray_r[0][0][0] = 1.0;
uarray_i[0][0][0] = 0.0;
for (int j = 1; j <= twojmax; j++) {
#if defined(_OPENMP)
#pragma omp for
#endif
for (int mb = 0; mb < j; mb++) {
uarray_r[j][0][mb] = 0.0;
uarray_i[j][0][mb] = 0.0;
for (int ma = 0; ma < j; ma++) {
rootpq = rootpqarray[j - ma][j - mb];
uarray_r[j][ma][mb] +=
rootpq *
(a_r * uarray_r[j - 1][ma][mb] +
a_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma][mb] +=
rootpq *
(a_r * uarray_i[j - 1][ma][mb] -
a_i * uarray_r[j - 1][ma][mb]);
rootpq = rootpqarray[ma + 1][j - mb];
uarray_r[j][ma + 1][mb] =
-rootpq *
(b_r * uarray_r[j - 1][ma][mb] +
b_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma + 1][mb] =
-rootpq *
(b_r * uarray_i[j - 1][ma][mb] -
b_i * uarray_r[j - 1][ma][mb]);
}
}
int mb = j;
uarray_r[j][0][mb] = 0.0;
uarray_i[j][0][mb] = 0.0;
#if defined(_OPENMP)
#pragma omp for
#endif
for (int ma = 0; ma < j; ma++) {
rootpq = rootpqarray[j - ma][mb];
uarray_r[j][ma][mb] +=
rootpq *
(b_r * uarray_r[j - 1][ma][mb - 1] -
b_i * uarray_i[j - 1][ma][mb - 1]);
uarray_i[j][ma][mb] +=
rootpq *
(b_r * uarray_i[j - 1][ma][mb - 1] +
b_i * uarray_r[j - 1][ma][mb - 1]);
rootpq = rootpqarray[ma + 1][mb];
uarray_r[j][ma + 1][mb] =
rootpq *
(a_r * uarray_r[j - 1][ma][mb - 1] -
a_i * uarray_i[j - 1][ma][mb - 1]);
uarray_i[j][ma + 1][mb] =
rootpq *
(a_r * uarray_i[j - 1][ma][mb - 1] +
a_i * uarray_r[j - 1][ma][mb - 1]);
}
}
}
/* ----------------------------------------------------------------------
compute derivatives of Wigner U-functions for one neighbor
see comments in compute_uarray()
------------------------------------------------------------------------- */
void SNA::compute_duarray(double x, double y, double z,
double z0, double r, double dz0dr,
double wj, double rcut)
{
double r0inv;
double a_r, a_i, b_r, b_i;
double da_r[3], da_i[3], db_r[3], db_i[3];
double dz0[3], dr0inv[3], dr0invdr;
double rootpq;
double rinv = 1.0 / r;
double ux = x * rinv;
double uy = y * rinv;
double uz = z * rinv;
r0inv = 1.0 / sqrt(r * r + z0 * z0);
a_r = z0 * r0inv;
a_i = -z * r0inv;
b_r = y * r0inv;
b_i = -x * r0inv;
dr0invdr = -pow(r0inv, 3.0) * (r + z0 * dz0dr);
dr0inv[0] = dr0invdr * ux;
dr0inv[1] = dr0invdr * uy;
dr0inv[2] = dr0invdr * uz;
dz0[0] = dz0dr * ux;
dz0[1] = dz0dr * uy;
dz0[2] = dz0dr * uz;
for (int k = 0; k < 3; k++) {
da_r[k] = dz0[k] * r0inv + z0 * dr0inv[k];
da_i[k] = -z * dr0inv[k];
}
da_i[2] += -r0inv;
for (int k = 0; k < 3; k++) {
db_r[k] = y * dr0inv[k];
db_i[k] = -x * dr0inv[k];
}
db_i[0] += -r0inv;
db_r[1] += r0inv;
uarray_r[0][0][0] = 1.0;
duarray_r[0][0][0][0] = 0.0;
duarray_r[0][0][0][1] = 0.0;
duarray_r[0][0][0][2] = 0.0;
uarray_i[0][0][0] = 0.0;
duarray_i[0][0][0][0] = 0.0;
duarray_i[0][0][0][1] = 0.0;
duarray_i[0][0][0][2] = 0.0;
for (int j = 1; j <= twojmax; j++) {
for (int mb = 0; 2*mb <= j; mb++) {
uarray_r[j][0][mb] = 0.0;
duarray_r[j][0][mb][0] = 0.0;
duarray_r[j][0][mb][1] = 0.0;
duarray_r[j][0][mb][2] = 0.0;
uarray_i[j][0][mb] = 0.0;
duarray_i[j][0][mb][0] = 0.0;
duarray_i[j][0][mb][1] = 0.0;
duarray_i[j][0][mb][2] = 0.0;
for (int ma = 0; ma < j; ma++) {
rootpq = rootpqarray[j - ma][j - mb];
uarray_r[j][ma][mb] += rootpq *
(a_r * uarray_r[j - 1][ma][mb] +
a_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma][mb] += rootpq *
(a_r * uarray_i[j - 1][ma][mb] -
a_i * uarray_r[j - 1][ma][mb]);
for (int k = 0; k < 3; k++) {
duarray_r[j][ma][mb][k] +=
rootpq * (da_r[k] * uarray_r[j - 1][ma][mb] +
da_i[k] * uarray_i[j - 1][ma][mb] +
a_r * duarray_r[j - 1][ma][mb][k] +
a_i * duarray_i[j - 1][ma][mb][k]);
duarray_i[j][ma][mb][k] +=
rootpq * (da_r[k] * uarray_i[j - 1][ma][mb] -
da_i[k] * uarray_r[j - 1][ma][mb] +
a_r * duarray_i[j - 1][ma][mb][k] -
a_i * duarray_r[j - 1][ma][mb][k]);
}
rootpq = rootpqarray[ma + 1][j - mb];
uarray_r[j][ma + 1][mb] =
-rootpq * (b_r * uarray_r[j - 1][ma][mb] +
b_i * uarray_i[j - 1][ma][mb]);
uarray_i[j][ma + 1][mb] =
-rootpq * (b_r * uarray_i[j - 1][ma][mb] -
b_i * uarray_r[j - 1][ma][mb]);
for (int k = 0; k < 3; k++) {
duarray_r[j][ma + 1][mb][k] =
-rootpq * (db_r[k] * uarray_r[j - 1][ma][mb] +
db_i[k] * uarray_i[j - 1][ma][mb] +
b_r * duarray_r[j - 1][ma][mb][k] +
b_i * duarray_i[j - 1][ma][mb][k]);
duarray_i[j][ma + 1][mb][k] =
-rootpq * (db_r[k] * uarray_i[j - 1][ma][mb] -
db_i[k] * uarray_r[j - 1][ma][mb] +
b_r * duarray_i[j - 1][ma][mb][k] -
b_i * duarray_r[j - 1][ma][mb][k]);
}
}
}
int mbpar = -1;
for (int mb = 0; 2*mb <= j; mb++) {
mbpar = -mbpar;
int mapar = -mbpar;
for (int ma = 0; ma <= j; ma++) {
mapar = -mapar;
if (mapar == 1) {
uarray_r[j][j-ma][j-mb] = uarray_r[j][ma][mb];
uarray_i[j][j-ma][j-mb] = -uarray_i[j][ma][mb];
for (int k = 0; k < 3; k++) {
duarray_r[j][j-ma][j-mb][k] = duarray_r[j][ma][mb][k];
duarray_i[j][j-ma][j-mb][k] = -duarray_i[j][ma][mb][k];
}
} else {
uarray_r[j][j-ma][j-mb] = -uarray_r[j][ma][mb];
uarray_i[j][j-ma][j-mb] = uarray_i[j][ma][mb];
for (int k = 0; k < 3; k++) {
duarray_r[j][j-ma][j-mb][k] = -duarray_r[j][ma][mb][k];
duarray_i[j][j-ma][j-mb][k] = duarray_i[j][ma][mb][k];
}
}
}
}
}
double sfac = compute_sfac(r, rcut);
double dsfac = compute_dsfac(r, rcut);
sfac *= wj;
dsfac *= wj;
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++)
for (int mb = 0; mb <= j; mb++) {
duarray_r[j][ma][mb][0] = dsfac * uarray_r[j][ma][mb] * ux +
sfac * duarray_r[j][ma][mb][0];
duarray_i[j][ma][mb][0] = dsfac * uarray_i[j][ma][mb] * ux +
sfac * duarray_i[j][ma][mb][0];
duarray_r[j][ma][mb][1] = dsfac * uarray_r[j][ma][mb] * uy +
sfac * duarray_r[j][ma][mb][1];
duarray_i[j][ma][mb][1] = dsfac * uarray_i[j][ma][mb] * uy +
sfac * duarray_i[j][ma][mb][1];
duarray_r[j][ma][mb][2] = dsfac * uarray_r[j][ma][mb] * uz +
sfac * duarray_r[j][ma][mb][2];
duarray_i[j][ma][mb][2] = dsfac * uarray_i[j][ma][mb] * uz +
sfac * duarray_i[j][ma][mb][2];
}
}
/* ----------------------------------------------------------------------
memory usage of arrays
------------------------------------------------------------------------- */
double SNA::memory_usage()
{
int jdim = twojmax + 1;
double bytes;
bytes = jdim * jdim * jdim * jdim * jdim * sizeof(double);
bytes += 2 * jdim * jdim * jdim * sizeof(complex<double>);
bytes += 2 * jdim * jdim * jdim * sizeof(double);
bytes += jdim * jdim * jdim * 3 * sizeof(complex<double>);
bytes += jdim * jdim * jdim * 3 * sizeof(double);
bytes += ncoeff * sizeof(double);
bytes += jdim * jdim * jdim * jdim * jdim * sizeof(complex<double>);
return bytes;
}
/* ---------------------------------------------------------------------- */
void SNA::create_twojmax_arrays()
{
int jdim = twojmax + 1;
memory->create(cgarray, jdim, jdim, jdim, jdim, jdim,
"sna:cgarray");
memory->create(rootpqarray, jdim+1, jdim+1,
"sna:rootpqarray");
memory->create(barray, jdim, jdim, jdim,
"sna:barray");
memory->create(dbarray, jdim, jdim, jdim, 3,
"sna:dbarray");
memory->create(duarray_r, jdim, jdim, jdim, 3,
"sna:duarray");
memory->create(duarray_i, jdim, jdim, jdim, 3,
"sna:duarray");
memory->create(uarray_r, jdim, jdim, jdim,
"sna:uarray");
memory->create(uarray_i, jdim, jdim, jdim,
"sna:uarray");
if(!use_shared_arrays) {
memory->create(uarraytot_r, jdim, jdim, jdim,
"sna:uarraytot");
memory->create(zarray_r, jdim, jdim, jdim, jdim, jdim,
"sna:zarray");
memory->create(uarraytot_i, jdim, jdim, jdim,
"sna:uarraytot");
memory->create(zarray_i, jdim, jdim, jdim, jdim, jdim,
"sna:zarray");
}
}
/* ---------------------------------------------------------------------- */
void SNA::destroy_twojmax_arrays()
{
memory->destroy(cgarray);
memory->destroy(rootpqarray);
memory->destroy(barray);
memory->destroy(dbarray);
memory->destroy(duarray_r);
memory->destroy(duarray_i);
memory->destroy(uarray_r);
memory->destroy(uarray_i);
if(!use_shared_arrays) {
memory->destroy(uarraytot_r);
memory->destroy(zarray_r);
memory->destroy(uarraytot_i);
memory->destroy(zarray_i);
}
}
/* ----------------------------------------------------------------------
factorial n, wrapper for precomputed table
------------------------------------------------------------------------- */
double SNA::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 SNA::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
};
/* ----------------------------------------------------------------------
the function delta given by VMK Eq. 8.2(1)
------------------------------------------------------------------------- */
double SNA::deltacg(int j1, int j2, int j)
{
double sfaccg = factorial((j1 + j2 + j) / 2 + 1);
return sqrt(factorial((j1 + j2 - j) / 2) *
factorial((j1 - j2 + j) / 2) *
factorial((-j1 + j2 + j) / 2) / sfaccg);
}
/* ----------------------------------------------------------------------
assign Clebsch-Gordan coefficients using
the quasi-binomial formula VMK 8.2.1(3)
------------------------------------------------------------------------- */
void SNA::init_clebsch_gordan()
{
double sum,dcg,sfaccg;
int m, aa2, bb2, cc2;
int ifac;
for (int j1 = 0; j1 <= twojmax; j1++)
for (int j2 = 0; j2 <= twojmax; j2++)
for (int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
for (int m1 = 0; m1 <= j1; m1 += 1) {
aa2 = 2 * m1 - j1;
for (int m2 = 0; m2 <= j2; m2 += 1) {
// -c <= cc <= c
bb2 = 2 * m2 - j2;
m = (aa2 + bb2 + j) / 2;
if(m < 0 || m > j) continue;
sum = 0.0;
for (int z = MAX(0, MAX(-(j - j2 + aa2)
/ 2, -(j - j1 - bb2) / 2));
z <= MIN((j1 + j2 - j) / 2,
MIN((j1 - aa2) / 2, (j2 + bb2) / 2));
z++) {
ifac = z % 2 ? -1 : 1;
sum += ifac /
(factorial(z) *
factorial((j1 + j2 - j) / 2 - z) *
factorial((j1 - aa2) / 2 - z) *
factorial((j2 + bb2) / 2 - z) *
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) *
factorial((j1 - aa2) / 2) *
factorial((j2 + bb2) / 2) *
factorial((j2 - bb2) / 2) *
factorial((j + cc2) / 2) *
factorial((j - cc2) / 2) *
(j + 1));
cgarray[j1][j2][j][m1][m2] = sum * dcg * sfaccg;
}
}
}
/* ----------------------------------------------------------------------
pre-compute table of sqrt[p/m2], p, q = 1,twojmax
the p = 0, q = 0 entries are allocated and skipped for convenience.
------------------------------------------------------------------------- */
void SNA::init_rootpqarray()
{
for (int p = 1; p <= twojmax; p++)
for (int q = 1; q <= twojmax; q++)
rootpqarray[p][q] = sqrt(static_cast<double>(p)/q);
}
/* ----------------------------------------------------------------------
a = j/2
------------------------------------------------------------------------- */
void SNA::jtostr(char* str, int j)
{
if(j % 2 == 0)
sprintf(str, "%d", j / 2);
else
sprintf(str, "%d/2", j);
}
/* ----------------------------------------------------------------------
aa = m - j/2
------------------------------------------------------------------------- */
void SNA::mtostr(char* str, int j, int m)
{
if(j % 2 == 0)
sprintf(str, "%d", m - j / 2);
else
sprintf(str, "%d/2", 2 * m - j);
}
/* ----------------------------------------------------------------------
list values of Clebsch-Gordan coefficients
using notation of VMK Table 8.11
------------------------------------------------------------------------- */
void SNA::print_clebsch_gordan(FILE* file)
{
char stra[20], strb[20], strc[20], straa[20], strbb[20], strcc[20];
int m, aa2, bb2;
fprintf(file, "a, aa, b, bb, c, cc, c(a,aa,b,bb,c,cc) \n");
for (int j1 = 0; j1 <= twojmax; j1++) {
jtostr(stra, j1);
for (int j2 = 0; j2 <= twojmax; j2++) {
jtostr(strb, j2);
for (int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
jtostr(strc, j);
for (int m1 = 0; m1 <= j1; m1 += 1) {
mtostr(straa, j1, m1);
aa2 = 2 * m1 - j1;
for (int m2 = 0; m2 <= j2; m2 += 1) {
bb2 = 2 * m2 - j2;
m = (aa2 + bb2 + j) / 2;
if(m < 0 || m > j) continue;
mtostr(strbb, j2, m2);
mtostr(strcc, j, m);
fprintf(file, "%s\t%s\t%s\t%s\t%s\t%s\t%g\n",
stra, straa, strb, strbb, strc, strcc,
cgarray[j1][j2][j][m1][m2]);
}
}
}
}
}
}
/* ---------------------------------------------------------------------- */
int SNA::compute_ncoeff()
{
int ncount;
ncount = 0;
for (int j1 = 0; j1 <= twojmax; j1++)
if(diagonalstyle == 0) {
for (int j2 = 0; j2 <= j1; j2++)
for (int j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2)
ncount++;
} else if(diagonalstyle == 1) {
int j2 = j1;
for (int j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2)
ncount++;
} else if(diagonalstyle == 2) {
ncount++;
} else if(diagonalstyle == 3) {
for (int j2 = 0; j2 <= j1; j2++)
for (int j = abs(j1 - j2);
j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) ncount++;
}
return ncount;
}
/* ---------------------------------------------------------------------- */
double SNA::compute_sfac(double r, double rcut)
{
if (switch_flag == 0) return 1.0;
if (switch_flag == 1) {
if(r <= rmin0) return 1.0;
else if(r > rcut) return 0.0;
else {
double rcutfac = MY_PI / (rcut - rmin0);
return 0.5 * (cos((r - rmin0) * rcutfac) + 1.0);
}
}
return 0.0;
}
/* ---------------------------------------------------------------------- */
double SNA::compute_dsfac(double r, double rcut)
{
if (switch_flag == 0) return 0.0;
if (switch_flag == 1) {
if(r <= rmin0) return 0.0;
else if(r > rcut) return 0.0;
else {
double rcutfac = MY_PI / (rcut - rmin0);
return -0.5 * sin((r - rmin0) * rcutfac) * rcutfac;
}
}
return 0.0;
}
Event Timeline
Log In to Comment