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

Reformat code with clang-format (Mozilla style guide)

Browse Source
This commit is contained in:
Sebastian Hütter
2017-06-09 18:05:17 +02:00
parent 7c7468ffc2
commit d3f31547f9
7 changed files with 2478 additions and 2211 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/USER-MEAMC/meam_cleanup.cpp
View File

@ -1,16 +1,17 @@
extern "C" {
#include "meam.h"
void meam_cleanup_(void) {
deallocate(meam_data.phirar6);
deallocate(meam_data.phirar5);
deallocate(meam_data.phirar4);
deallocate(meam_data.phirar3);
deallocate(meam_data.phirar2);
deallocate(meam_data.phirar1);
deallocate(meam_data.phirar);
deallocate(meam_data.phir);
void
meam_cleanup_(void)
{
deallocate(meam_data.phirar6);
deallocate(meam_data.phirar5);
deallocate(meam_data.phirar4);
deallocate(meam_data.phirar3);
deallocate(meam_data.phirar2);
deallocate(meam_data.phirar1);
deallocate(meam_data.phirar);
deallocate(meam_data.phir);
}
}

467
src/USER-MEAMC/meam_dens_final.cpp
View File

@ -1,12 +1,14 @@
extern "C"{
extern "C" {
#include "meam.h"
#include <math.h>
void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag);
void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG);
void G_gam(double Gamma, int ibar, double gsmooth_factor, double* G,
int* errorflag);
void dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G,
double* dG);
// in meam_setup_done
void get_shpfcn(double *s /* s(3) */, lattice_t latt);
void get_shpfcn(double* s /* s(3) */, lattice_t latt);
// Extern "C" declaration has the form:
//
@ -25,241 +27,266 @@ extern "C"{
// dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
//
void meam_dens_final_(int *nlocal, int *nmax,
int *eflag_either, int *eflag_global, int *eflag_atom, double *eng_vdwl, double *eatom,
int *ntype, int *type, int *fmap,
double *Arho1, double *Arho2, double *Arho2b, double *Arho3, double *Arho3b, double *t_ave, double *tsq_ave,
double *Gamma, double *dGamma1, double *dGamma2, double *dGamma3,
double *rho, double *rho0, double *rho1, double *rho2, double *rho3, double *fp, int *errorflag)
{
arrdim2v(Arho1,3,*nmax)
arrdim2v(Arho2,6,*nmax)
arrdim2v(Arho3,10,*nmax)
arrdim2v(Arho3b,3,*nmax)
arrdim2v(t_ave,3,*nmax)
arrdim2v(tsq_ave,3,*nmax)
void
meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype,
int* type, int* fmap, double* Arho1, double* Arho2,
double* Arho2b, double* Arho3, double* Arho3b, double* t_ave,
double* tsq_ave, double* Gamma, double* dGamma1,
double* dGamma2, double* dGamma3, double* rho, double* rho0,
double* rho1, double* rho2, double* rho3, double* fp,
int* errorflag)
{
int i, elti;
int m;
double rhob, G, dG, Gbar, dGbar, gam, shp[3+1], Z;
double B, denom, rho_bkgd;
arrdim2v(Arho1, 3, *nmax);
arrdim2v(Arho2, 6, *nmax);
arrdim2v(Arho3, 10, *nmax);
arrdim2v(Arho3b, 3, *nmax);
arrdim2v(t_ave, 3, *nmax);
arrdim2v(tsq_ave, 3, *nmax);
// Complete the calculation of density
int i, elti;
int m;
double rhob, G, dG, Gbar, dGbar, gam, shp[3 + 1], Z;
double B, denom, rho_bkgd;
for(i=1; i<=*nlocal; i++) {
elti = arr1v(fmap,arr1v(type,i));
if (elti > 0) {
arr1v(rho1,i) = 0.0;
arr1v(rho2,i) = -1.0/3.0*arr1v(Arho2b,i)*arr1v(Arho2b,i);
arr1v(rho3,i) = 0.0;
for(m=1; m<=3; m++) {
arr1v(rho1,i) = arr1v(rho1,i) + arr2v(Arho1,m,i)*arr2v(Arho1,m,i);
arr1v(rho3,i) = arr1v(rho3,i) - 3.0/5.0*arr2v(Arho3b,m,i)*arr2v(Arho3b,m,i);
}
for(m=1; m<=6; m++) {
arr1v(rho2,i) = arr1v(rho2,i) + meam_data.v2D[m]*arr2v(Arho2,m,i)*arr2v(Arho2,m,i);
}
for(m=1; m<=10; m++) {
arr1v(rho3,i) = arr1v(rho3,i) + meam_data.v3D[m]*arr2v(Arho3,m,i)*arr2v(Arho3,m,i);
}
// Complete the calculation of density
if( arr1v(rho0,i) > 0.0 ) {
if (meam_data.ialloy==1) {
arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr2v(tsq_ave,1,i);
arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr2v(tsq_ave,2,i);
arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr2v(tsq_ave,3,i);
} else if (meam_data.ialloy==2) {
arr2v(t_ave,1,i) = meam_data.t1_meam[elti];
arr2v(t_ave,2,i) = meam_data.t2_meam[elti];
arr2v(t_ave,3,i) = meam_data.t3_meam[elti];
} else {
arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr1v(rho0,i);
arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr1v(rho0,i);
arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr1v(rho0,i);
}
}
for (i = 1; i <= *nlocal; i++) {
elti = arr1v(fmap, arr1v(type, i));
if (elti > 0) {
arr1v(rho1, i) = 0.0;
arr1v(rho2, i) = -1.0 / 3.0 * arr1v(Arho2b, i) * arr1v(Arho2b, i);
arr1v(rho3, i) = 0.0;
for (m = 1; m <= 3; m++) {
arr1v(rho1, i) =
arr1v(rho1, i) + arr2v(Arho1, m, i) * arr2v(Arho1, m, i);
arr1v(rho3, i) = arr1v(rho3, i) -
3.0 / 5.0 * arr2v(Arho3b, m, i) * arr2v(Arho3b, m, i);
}
for (m = 1; m <= 6; m++) {
arr1v(rho2, i) =
arr1v(rho2, i) +
meam_data.v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i);
}
for (m = 1; m <= 10; m++) {
arr1v(rho3, i) =
arr1v(rho3, i) +
meam_data.v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i);
}
arr1v(Gamma,i) = arr2v(t_ave,1,i)*arr1v(rho1,i) + arr2v(t_ave,2,i)*arr1v(rho2,i) + arr2v(t_ave,3,i)*arr1v(rho3,i);
if( arr1v(rho0,i) > 0.0 ) {
arr1v(Gamma,i) = arr1v(Gamma,i)/(arr1v(rho0,i)*arr1v(rho0,i));
}
Z = meam_data.Z_meam[elti];
G_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti], meam_data.gsmooth_factor,&G,errorflag);
if (*errorflag!=0) return;
get_shpfcn(shp,meam_data.lattce_meam[elti][elti]);
if (meam_data.ibar_meam[elti]<=0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
if (meam_data.mix_ref_t==1) {
gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z);
} else {
gam = (meam_data.t1_meam[elti]*shp[1]+meam_data.t2_meam[elti]*shp[2] +meam_data.t3_meam[elti]*shp[3])/(Z*Z);
}
G_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&Gbar,errorflag);
}
arr1v(rho,i) = arr1v(rho0,i) * G;
if (meam_data.mix_ref_t==1) {
if (meam_data.ibar_meam[elti]<=0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z);
dG_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor, &Gbar,&dGbar);
}
rho_bkgd = meam_data.rho0_meam[elti]*Z*Gbar;
} else {
if (meam_data.bkgd_dyn==1) {
rho_bkgd = meam_data.rho0_meam[elti]*Z;
} else {
rho_bkgd = meam_data.rho_ref_meam[elti];
}
}
rhob = arr1v(rho,i)/rho_bkgd;
denom = 1.0/rho_bkgd;
dG_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&G,&dG);
arr1v(dGamma1,i) = (G - 2*dG*arr1v(Gamma,i))*denom;
if( !iszero(arr1v(rho0,i)) ) {
arr1v(dGamma2,i) = (dG/arr1v(rho0,i))*denom;
} else {
arr1v(dGamma2,i) = 0.0;
}
// dGamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (meam_data.mix_ref_t==1) {
arr1v(dGamma3,i) = arr1v(rho0,i)*G*dGbar/(Gbar*Z*Z)*denom;
} else {
arr1v(dGamma3,i) = 0.0;
}
B = meam_data.A_meam[elti]*meam_data.Ec_meam[elti][elti];
if( !iszero(rhob) ) {
if (meam_data.emb_lin_neg==1 && rhob<=0) {
arr1v(fp,i) = -B;
} else {
arr1v(fp,i) = B*(log(rhob)+1.0);
}
if (*eflag_either!=0) {
if (*eflag_global!=0) {
if (meam_data.emb_lin_neg==1 && rhob<=0) {
*eng_vdwl = *eng_vdwl - B*rhob;
} else {
*eng_vdwl = *eng_vdwl + B*rhob*log(rhob);
}
}
if (*eflag_atom!=0) {
if (meam_data.emb_lin_neg==1 && rhob<=0) {
arr1v(eatom,i) = arr1v(eatom,i) - B*rhob;
} else {
arr1v(eatom,i) = arr1v(eatom,i) + B*rhob*log(rhob);
}
}
}
} else {
if (meam_data.emb_lin_neg==1) {
arr1v(fp,i) = -B;
} else {
arr1v(fp,i) = B;
}
}
if (arr1v(rho0, i) > 0.0) {
if (meam_data.ialloy == 1) {
arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr2v(tsq_ave, 1, i);
arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr2v(tsq_ave, 2, i);
arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr2v(tsq_ave, 3, i);
} else if (meam_data.ialloy == 2) {
arr2v(t_ave, 1, i) = meam_data.t1_meam[elti];
arr2v(t_ave, 2, i) = meam_data.t2_meam[elti];
arr2v(t_ave, 3, i) = meam_data.t3_meam[elti];
} else {
arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr1v(rho0, i);
arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr1v(rho0, i);
arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr1v(rho0, i);
}
}
arr1v(Gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) +
arr2v(t_ave, 2, i) * arr1v(rho2, i) +
arr2v(t_ave, 3, i) * arr1v(rho3, i);
if (arr1v(rho0, i) > 0.0) {
arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i));
}
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
Z = meam_data.Z_meam[elti];
void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag)
{
// Compute G(Gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = exp(Gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
double gsmooth_switchpoint;
if (ibar==0 || ibar==4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1);
if (Gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
*G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma),gsmooth_factor);
*G = sqrt(*G);
G_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti],
meam_data.gsmooth_factor, &G, errorflag);
if (*errorflag != 0)
return;
get_shpfcn(shp, meam_data.lattce_meam[elti][elti]);
if (meam_data.ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
if (meam_data.mix_ref_t == 1) {
gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] +
arr2v(t_ave, 3, i) * shp[3]) /
(Z * Z);
} else {
*G = sqrt(1.0+Gamma);
gam = (meam_data.t1_meam[elti] * shp[1] +
meam_data.t2_meam[elti] * shp[2] +
meam_data.t3_meam[elti] * shp[3]) /
(Z * Z);
}
} else if (ibar==1) {
*G = fm_exp(Gamma/2.0);
} else if (ibar==3) {
*G = 2.0/(1.0+exp(-Gamma));
} else if (ibar==-5) {
if ((1.0+Gamma)>=0) {
*G = sqrt(1.0+Gamma);
G_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor, &Gbar,
errorflag);
}
arr1v(rho, i) = arr1v(rho0, i) * G;
if (meam_data.mix_ref_t == 1) {
if (meam_data.ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
*G = -sqrt(-1.0-Gamma);
gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] +
arr2v(t_ave, 3, i) * shp[3]) /
(Z * Z);
dG_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor,
&Gbar, &dGbar);
}
rho_bkgd = meam_data.rho0_meam[elti] * Z * Gbar;
} else {
if (meam_data.bkgd_dyn == 1) {
rho_bkgd = meam_data.rho0_meam[elti] * Z;
} else {
rho_bkgd = meam_data.rho_ref_meam[elti];
}
}
rhob = arr1v(rho, i) / rho_bkgd;
denom = 1.0 / rho_bkgd;
dG_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti],
meam_data.gsmooth_factor, &G, &dG);
arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom;
if (!iszero(arr1v(rho0, i))) {
arr1v(dGamma2, i) = (dG / arr1v(rho0, i)) * denom;
} else {
arr1v(dGamma2, i) = 0.0;
}
// dGamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (meam_data.mix_ref_t == 1) {
arr1v(dGamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom;
} else {
arr1v(dGamma3, i) = 0.0;
}
B = meam_data.A_meam[elti] * meam_data.Ec_meam[elti][elti];
if (!iszero(rhob)) {
if (meam_data.emb_lin_neg == 1 && rhob <= 0) {
arr1v(fp, i) = -B;
} else {
arr1v(fp, i) = B * (log(rhob) + 1.0);
}
if (*eflag_either != 0) {
if (*eflag_global != 0) {
if (meam_data.emb_lin_neg == 1 && rhob <= 0) {
*eng_vdwl = *eng_vdwl - B * rhob;
} else {
*eng_vdwl = *eng_vdwl + B * rhob * log(rhob);
}
}
if (*eflag_atom != 0) {
if (meam_data.emb_lin_neg == 1 && rhob <= 0) {
arr1v(eatom, i) = arr1v(eatom, i) - B * rhob;
} else {
arr1v(eatom, i) = arr1v(eatom, i) + B * rhob * log(rhob);
}
}
}
} else {
*errorflag = 1;
}
}
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG)
{
// Compute G(Gamma) and dG(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = fm_exp(Gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+fm_exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
double gsmooth_switchpoint;
if (ibar==0 || ibar==4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1);
if (Gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
*G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma), gsmooth_factor);
*G = sqrt(*G);
*dG = -gsmooth_factor* *G/(2.0*Gamma);
if (meam_data.emb_lin_neg == 1) {
arr1v(fp, i) = -B;
} else {
*G = sqrt(1.0+Gamma);
*dG = 1.0/(2.0* *G);
}
} else if (ibar==1) {
*G = fm_exp(Gamma/2.0);
*dG = *G/2.0;
} else if (ibar==3) {
*G = 2.0/(1.0+fm_exp(-Gamma));
*dG = *G*(2.0-*G)/2;
} else if (ibar==-5) {
if ((1.0+Gamma)>=0) {
*G = sqrt(1.0+Gamma);
*dG = 1.0/(2.0* *G);
} else {
*G = -sqrt(-1.0-Gamma);
*dG = -1.0/(2.0* *G);
arr1v(fp, i) = B;
}
}
}
}
}
}
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag)
{
// Compute G(Gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = exp(Gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
double gsmooth_switchpoint;
if (ibar == 0 || ibar == 4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (Gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / Gamma), gsmooth_factor);
*G = sqrt(*G);
} else {
*G = sqrt(1.0 + Gamma);
}
} else if (ibar == 1) {
*G = fm_exp(Gamma / 2.0);
} else if (ibar == 3) {
*G = 2.0 / (1.0 + exp(-Gamma));
} else if (ibar == -5) {
if ((1.0 + Gamma) >= 0) {
*G = sqrt(1.0 + Gamma);
} else {
*G = -sqrt(-1.0 - Gamma);
}
} else {
*errorflag = 1;
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG)
{
// Compute G(Gamma) and dG(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = fm_exp(Gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+fm_exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
double gsmooth_switchpoint;
if (ibar == 0 || ibar == 4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (Gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / Gamma), gsmooth_factor);
*G = sqrt(*G);
*dG = -gsmooth_factor * *G / (2.0 * Gamma);
} else {
*G = sqrt(1.0 + Gamma);
*dG = 1.0 / (2.0 * *G);
}
} else if (ibar == 1) {
*G = fm_exp(Gamma / 2.0);
*dG = *G / 2.0;
} else if (ibar == 3) {
*G = 2.0 / (1.0 + fm_exp(-Gamma));
*dG = *G * (2.0 - *G) / 2;
} else if (ibar == -5) {
if ((1.0 + Gamma) >= 0) {
*G = sqrt(1.0 + Gamma);
*dG = 1.0 / (2.0 * *G);
} else {
*G = -sqrt(-1.0 - Gamma);
*dG = -1.0 / (2.0 * *G);
}
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
}

890
src/USER-MEAMC/meam_dens_init.cpp
View File

@ -2,25 +2,30 @@ extern "C" {
#include "meam.h"
#include <math.h>
void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x,
int numneigh, int *firstneigh,
int numneigh_full, int *firstneigh_full,
int ntype, int *type, int *fmap);
void calc_rho1(int i, int nmax, int ntype, int *type, int *fmap, double *x,
int numneigh, int *firstneigh,
double *scrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
double *arho3, double *arho3b, double *t_ave, double *tsq_ave);
void screen(int i, int j, int nmax, double *x, double rijsq, double *sij, int numneigh_full, int *firstneigh_full, int ntype, int *type, int *fmap);
void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair);
void fcut(double xi, double * fc);
void dfcut(double xi, double *fc, double *dfc);
void dCfunc(double rij2,double rik2,double rjk2,double *dCikj);
void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2);
void getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair,
double* x, int numneigh, int* firstneigh, int numneigh_full,
int* firstneigh_full, int ntype, int* type, int* fmap);
void calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair,
double* rho0, double* arho1, double* arho2, double* arho2b,
double* arho3, double* arho3b, double* t_ave, double* tsq_ave);
void screen(int i, int j, int nmax, double* x, double rijsq, double* sij,
int numneigh_full, int* firstneigh_full, int ntype, int* type,
int* fmap);
void dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2,
double* dsij1, double* dsij2, int ntype, int* type, int* fmap,
double* x, double* scrfcn, double* fcpair);
void fcut(double xi, double* fc);
void dfcut(double xi, double* fc, double* dfc);
void dCfunc(double rij2, double rik2, double rjk2, double* dCikj);
void dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1,
double* dCikj2);
// Extern "C" declaration has the form:
//
// void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *,
// void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *,
// double *,
// int *, int *, int *, int *,
// double *, double *, double *, double *, double *, double *,
// double *, double *, double *, double *, double *, int *);
@ -35,470 +40,505 @@ extern "C" {
// &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
//
void meam_dens_init_(int *i, int *nmax,
int *ntype, int *type, int *fmap, double *x,
int *numneigh, int *firstneigh,
int *numneigh_full, int *firstneigh_full,
double *scrfcn, double *dscrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
double *arho3, double *arho3b, double *t_ave, double *tsq_ave, int *errorflag)
{
*errorflag = 0;
void
meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x,
int* numneigh, int* firstneigh, int* numneigh_full,
int* firstneigh_full, double* scrfcn, double* dscrfcn,
double* fcpair, double* rho0, double* arho1, double* arho2,
double* arho2b, double* arho3, double* arho3b, double* t_ave,
double* tsq_ave, int* errorflag)
{
*errorflag = 0;
// Compute screening function and derivatives
getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x,
*numneigh, firstneigh,
*numneigh_full, firstneigh_full,
*ntype, type, fmap);
// Compute screening function and derivatives
getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x, *numneigh, firstneigh,
*numneigh_full, firstneigh_full, *ntype, type, fmap);
// Calculate intermediate density terms to be communicated
calc_rho1(*i, *nmax, *ntype, type, fmap, x,
*numneigh, firstneigh,
scrfcn, fcpair, rho0, arho1, arho2, arho2b,
arho3, arho3b, t_ave, tsq_ave);
// Calculate intermediate density terms to be communicated
calc_rho1(*i, *nmax, *ntype, type, fmap, x, *numneigh, firstneigh, scrfcn,
fcpair, rho0, arho1, arho2, arho2b, arho3, arho3b, t_ave, tsq_ave);
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair,
double* x, int numneigh, int* firstneigh, int numneigh_full,
int* firstneigh_full, int ntype, int* type, int* fmap)
{
void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x,
int numneigh, int *firstneigh,
int numneigh_full, int *firstneigh_full,
int ntype, int *type, int *fmap)
{
arrdim2v(x,3,nmax)
arrdim2v(x, 3, nmax);
int jn,j,kn,k;
int elti,eltj,eltk;
double xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij;
double xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2/*,rik*/;
double xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2/*,rjk*/;
double xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj;
double Cmin,Cmax,delc,/*ebound,*/rbound,a,coef1,coef2;
//double coef1a,coef1b,coef2a,coef2b;
double dCikj;
/*unused:double dC1a,dC1b,dC2a,dC2b;*/
double rnorm,fc,dfc,drinv;
int jn, j, kn, k;
int elti, eltj, eltk;
double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij;
double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/;
double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/;
double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj;
double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2;
// double coef1a,coef1b,coef2a,coef2b;
double dCikj;
/*unused:double dC1a,dC1b,dC2a,dC2b;*/
double rnorm, fc, dfc, drinv;
drinv = 1.0/meam_data.delr_meam;
elti = arr1v(fmap,arr1v(type,i));
drinv = 1.0 / meam_data.delr_meam;
elti = arr1v(fmap, arr1v(type, i));
if (elti > 0) {
if (elti > 0) {
xitmp = arr2v(x,1,i);
yitmp = arr2v(x,2,i);
zitmp = arr2v(x,3,i);
xitmp = arr2v(x, 1, i);
yitmp = arr2v(x, 2, i);
zitmp = arr2v(x, 3, i);
for(jn=1; jn<=numneigh; jn++) {
j = arr1v(firstneigh,jn);
for (jn = 1; jn <= numneigh; jn++) {
j = arr1v(firstneigh, jn);
eltj = arr1v(fmap,arr1v(type,j));
if (eltj > 0) {
// First compute screening function itself, sij
xjtmp = arr2v(x,1,j);
yjtmp = arr2v(x,2,j);
zjtmp = arr2v(x,3,j);
delxij = xjtmp - xitmp;
delyij = yjtmp - yitmp;
delzij = zjtmp - zitmp;
rij2 = delxij*delxij + delyij*delyij + delzij*delzij;
rij = sqrt(rij2);
if (rij > meam_data.rc_meam) {
fcij = 0.0;
dfcij = 0.0;
sij = 0.0;
} else {
rnorm = (meam_data.rc_meam-rij)*drinv;
screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap);
dfcut(rnorm,&fc,&dfc);
fcij = fc;
dfcij = dfc*drinv;
}
// Now compute derivatives
arr1v(dscrfcn,jn) = 0.0;
sfcij = sij*fcij;
if (iszero(sfcij) || iszero(sfcij-1.0)) goto LABEL_100;
rbound = meam_data.ebound_meam[elti][eltj] * rij2;
for(kn=1; kn<=numneigh_full; kn++) {
k = arr1v(firstneigh_full,kn);
if (k==j) continue;
eltk = arr1v(fmap,arr1v(type,k));
if (eltk==0) continue;
xktmp = arr2v(x,1,k);
yktmp = arr2v(x,2,k);
zktmp = arr2v(x,3,k);
delxjk = xktmp - xjtmp;
delyjk = yktmp - yjtmp;
delzjk = zktmp - zjtmp;
rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk;
if (rjk2 > rbound) continue;
delxik = xktmp - xitmp;
delyik = yktmp - yitmp;
delzik = zktmp - zitmp;
rik2 = delxik*delxik + delyik*delyik + delzik*delzik;
if (rik2 > rbound) continue;
xik = rik2/rij2;
xjk = rjk2/rij2;
a = 1 - (xik-xjk)*(xik-xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a<=0.0) continue;
cikj = (2.0*(xik+xjk) + a - 2.0)/a;
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
if (cikj>=Cmax) {
continue;
// Note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case here
// (other negative cikj cases were handled by the test on "a" above)
// Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
} else {
delc = Cmax - Cmin;
cikj = (cikj-Cmin)/delc;
dfcut(cikj,&sikj,&dfikj);
coef1 = dfikj/(delc*sikj);
dCfunc(rij2,rik2,rjk2,&dCikj);
arr1v(dscrfcn,jn) = arr1v(dscrfcn,jn) + coef1*dCikj;
}
}
coef1 = sfcij;
coef2 = sij*dfcij/rij;
arr1v(dscrfcn,jn) = arr1v(dscrfcn,jn)*coef1 - coef2;
LABEL_100:
arr1v(scrfcn,jn) = sij;
arr1v(fcpair,jn) = fcij;
}
eltj = arr1v(fmap, arr1v(type, j));
if (eltj > 0) {
// First compute screening function itself, sij
xjtmp = arr2v(x, 1, j);
yjtmp = arr2v(x, 2, j);
zjtmp = arr2v(x, 3, j);
delxij = xjtmp - xitmp;
delyij = yjtmp - yitmp;
delzij = zjtmp - zitmp;
rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
rij = sqrt(rij2);
if (rij > meam_data.rc_meam) {
fcij = 0.0;
dfcij = 0.0;
sij = 0.0;
} else {
rnorm = (meam_data.rc_meam - rij) * drinv;
screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full,
ntype, type, fmap);
dfcut(rnorm, &fc, &dfc);
fcij = fc;
dfcij = dfc * drinv;
}
// Now compute derivatives
arr1v(dscrfcn, jn) = 0.0;
sfcij = sij * fcij;
if (iszero(sfcij) || iszero(sfcij - 1.0))
goto LABEL_100;
rbound = meam_data.ebound_meam[elti][eltj] * rij2;
for (kn = 1; kn <= numneigh_full; kn++) {
k = arr1v(firstneigh_full, kn);
if (k == j)
continue;
eltk = arr1v(fmap, arr1v(type, k));
if (eltk == 0)
continue;
xktmp = arr2v(x, 1, k);
yktmp = arr2v(x, 2, k);
zktmp = arr2v(x, 3, k);
delxjk = xktmp - xjtmp;
delyjk = yktmp - yjtmp;
delzjk = zktmp - zjtmp;
rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
if (rjk2 > rbound)
continue;
delxik = xktmp - xitmp;
delyik = yktmp - yitmp;
delzik = zktmp - zitmp;
rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
if (rik2 > rbound)
continue;
xik = rik2 / rij2;
xjk = rjk2 / rij2;
a = 1 - (xik - xjk) * (xik - xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a <= 0.0)
continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax) {
continue;
// Note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case
// here
// (other negative cikj cases were handled by the test on "a"
// above)
// Note that we never have 0<cikj<Cmin here, else sij=0
// (rejected above)
} else {
delc = Cmax - Cmin;
cikj = (cikj - Cmin) / delc;
dfcut(cikj, &sikj, &dfikj);
coef1 = dfikj / (delc * sikj);
dCfunc(rij2, rik2, rjk2, &dCikj);
arr1v(dscrfcn, jn) = arr1v(dscrfcn, jn) + coef1 * dCikj;
}
}
coef1 = sfcij;
coef2 = sij * dfcij / rij;
arr1v(dscrfcn, jn) = arr1v(dscrfcn, jn) * coef1 - coef2;
LABEL_100:
arr1v(scrfcn, jn) = sij;
arr1v(fcpair, jn) = fcij;
}
}
}
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair,
double* rho0, double* arho1, double* arho2, double* arho2b,
double* arho3, double* arho3b, double* t_ave, double* tsq_ave)
{
arrdim2v(x, 3, nmax);
arrdim2v(arho1, 3, nmax);
arrdim2v(arho2, 6, nmax);
arrdim2v(arho3, 10, nmax);
arrdim2v(arho3b, 3, nmax);
arrdim2v(t_ave, 3, nmax);
arrdim2v(tsq_ave, 3, nmax);
void calc_rho1(int i, int nmax, int ntype, int *type, int *fmap, double *x,
int numneigh, int *firstneigh,
double *scrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
double *arho3, double *arho3b, double *t_ave, double *tsq_ave)
{
arrdim2v(x,3,nmax)
arrdim2v(arho1,3,nmax)
arrdim2v(arho2,6,nmax)
arrdim2v(arho3,10,nmax)
arrdim2v(arho3b,3,nmax)
arrdim2v(t_ave,3,nmax)
arrdim2v(tsq_ave,3,nmax)
int jn, j, m, n, p, elti, eltj;
int nv2, nv3;
double xtmp, ytmp, ztmp, delij[3 + 1], rij2, rij, sij;
double ai, aj, rhoa0j, rhoa1j, rhoa2j, rhoa3j, A1j, A2j, A3j;
// double G,Gbar,gam,shp[3+1];
double ro0i, ro0j;
double rhoa0i, rhoa1i, rhoa2i, rhoa3i, A1i, A2i, A3i;
int jn,j,m,n,p,elti,eltj;
int nv2,nv3;
double xtmp,ytmp,ztmp,delij[3+1],rij2,rij,sij;
double ai,aj,rhoa0j,rhoa1j,rhoa2j,rhoa3j,A1j,A2j,A3j;
//double G,Gbar,gam,shp[3+1];
double ro0i,ro0j;
double rhoa0i,rhoa1i,rhoa2i,rhoa3i,A1i,A2i,A3i;
elti = arr1v(fmap, arr1v(type, i));
xtmp = arr2v(x, 1, i);
ytmp = arr2v(x, 2, i);
ztmp = arr2v(x, 3, i);
for (jn = 1; jn <= numneigh; jn++) {
if (!iszero(arr1v(scrfcn, jn))) {
j = arr1v(firstneigh, jn);
sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn);
delij[1] = arr2v(x, 1, j) - xtmp;
delij[2] = arr2v(x, 2, j) - ytmp;
delij[3] = arr2v(x, 3, j) - ztmp;
rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3];
if (rij2 < meam_data.cutforcesq) {
eltj = arr1v(fmap, arr1v(type, j));
rij = sqrt(rij2);
ai = rij / meam_data.re_meam[elti][elti] - 1.0;
aj = rij / meam_data.re_meam[eltj][eltj] - 1.0;
ro0i = meam_data.rho0_meam[elti];
ro0j = meam_data.rho0_meam[eltj];
rhoa0j = ro0j * fm_exp(-meam_data.beta0_meam[eltj] * aj) * sij;
rhoa1j = ro0j * fm_exp(-meam_data.beta1_meam[eltj] * aj) * sij;
rhoa2j = ro0j * fm_exp(-meam_data.beta2_meam[eltj] * aj) * sij;
rhoa3j = ro0j * fm_exp(-meam_data.beta3_meam[eltj] * aj) * sij;
rhoa0i = ro0i * fm_exp(-meam_data.beta0_meam[elti] * ai) * sij;
rhoa1i = ro0i * fm_exp(-meam_data.beta1_meam[elti] * ai) * sij;
rhoa2i = ro0i * fm_exp(-meam_data.beta2_meam[elti] * ai) * sij;
rhoa3i = ro0i * fm_exp(-meam_data.beta3_meam[elti] * ai) * sij;
if (meam_data.ialloy == 1) {
rhoa1j = rhoa1j * meam_data.t1_meam[eltj];
rhoa2j = rhoa2j * meam_data.t2_meam[eltj];
rhoa3j = rhoa3j * meam_data.t3_meam[eltj];
rhoa1i = rhoa1i * meam_data.t1_meam[elti];
rhoa2i = rhoa2i * meam_data.t2_meam[elti];
rhoa3i = rhoa3i * meam_data.t3_meam[elti];
}
arr1v(rho0, i) = arr1v(rho0, i) + rhoa0j;
arr1v(rho0, j) = arr1v(rho0, j) + rhoa0i;
// For ialloy = 2, use single-element value (not average)
if (meam_data.ialloy != 2) {
arr2v(t_ave, 1, i) =
arr2v(t_ave, 1, i) + meam_data.t1_meam[eltj] * rhoa0j;
arr2v(t_ave, 2, i) =
arr2v(t_ave, 2, i) + meam_data.t2_meam[eltj] * rhoa0j;
arr2v(t_ave, 3, i) =
arr2v(t_ave, 3, i) + meam_data.t3_meam[eltj] * rhoa0j;
arr2v(t_ave, 1, j) =
arr2v(t_ave, 1, j) + meam_data.t1_meam[elti] * rhoa0i;
arr2v(t_ave, 2, j) =
arr2v(t_ave, 2, j) + meam_data.t2_meam[elti] * rhoa0i;
arr2v(t_ave, 3, j) =
arr2v(t_ave, 3, j) + meam_data.t3_meam[elti] * rhoa0i;
}
if (meam_data.ialloy == 1) {
arr2v(tsq_ave, 1, i) =
arr2v(tsq_ave, 1, i) +
meam_data.t1_meam[eltj] * meam_data.t1_meam[eltj] * rhoa0j;
arr2v(tsq_ave, 2, i) =
arr2v(tsq_ave, 2, i) +
meam_data.t2_meam[eltj] * meam_data.t2_meam[eltj] * rhoa0j;
arr2v(tsq_ave, 3, i) =
arr2v(tsq_ave, 3, i) +
meam_data.t3_meam[eltj] * meam_data.t3_meam[eltj] * rhoa0j;
arr2v(tsq_ave, 1, j) =
arr2v(tsq_ave, 1, j) +
meam_data.t1_meam[elti] * meam_data.t1_meam[elti] * rhoa0i;
arr2v(tsq_ave, 2, j) =
arr2v(tsq_ave, 2, j) +
meam_data.t2_meam[elti] * meam_data.t2_meam[elti] * rhoa0i;
arr2v(tsq_ave, 3, j) =
arr2v(tsq_ave, 3, j) +
meam_data.t3_meam[elti] * meam_data.t3_meam[elti] * rhoa0i;
}
arr1v(arho2b, i) = arr1v(arho2b, i) + rhoa2j;
arr1v(arho2b, j) = arr1v(arho2b, j) + rhoa2i;
elti = arr1v(fmap,arr1v(type,i));
xtmp = arr2v(x,1,i);
ytmp = arr2v(x,2,i);
ztmp = arr2v(x,3,i);
for(jn=1; jn<=numneigh; jn++) {
if (!iszero(arr1v(scrfcn,jn))) {
j = arr1v(firstneigh,jn);
sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn);
delij[1] = arr2v(x,1,j) - xtmp;
delij[2] = arr2v(x,2,j) - ytmp;
delij[3] = arr2v(x,3,j) - ztmp;
rij2 = delij[1]*delij[1] + delij[2]*delij[2] + delij[3]*delij[3];
if (rij2 < meam_data.cutforcesq) {
eltj = arr1v(fmap,arr1v(type,j));
rij = sqrt(rij2);
ai = rij/meam_data.re_meam[elti][elti] - 1.0;
aj = rij/meam_data.re_meam[eltj][eltj] - 1.0;
ro0i = meam_data.rho0_meam[elti];
ro0j = meam_data.rho0_meam[eltj];
rhoa0j = ro0j*fm_exp(-meam_data.beta0_meam[eltj]*aj)*sij;
rhoa1j = ro0j*fm_exp(-meam_data.beta1_meam[eltj]*aj)*sij;
rhoa2j = ro0j*fm_exp(-meam_data.beta2_meam[eltj]*aj)*sij;
rhoa3j = ro0j*fm_exp(-meam_data.beta3_meam[eltj]*aj)*sij;
rhoa0i = ro0i*fm_exp(-meam_data.beta0_meam[elti]*ai)*sij;
rhoa1i = ro0i*fm_exp(-meam_data.beta1_meam[elti]*ai)*sij;
rhoa2i = ro0i*fm_exp(-meam_data.beta2_meam[elti]*ai)*sij;
rhoa3i = ro0i*fm_exp(-meam_data.beta3_meam[elti]*ai)*sij;
if (meam_data.ialloy==1) {
rhoa1j = rhoa1j * meam_data.t1_meam[eltj];
rhoa2j = rhoa2j * meam_data.t2_meam[eltj];
rhoa3j = rhoa3j * meam_data.t3_meam[eltj];
rhoa1i = rhoa1i * meam_data.t1_meam[elti];
rhoa2i = rhoa2i * meam_data.t2_meam[elti];
rhoa3i = rhoa3i * meam_data.t3_meam[elti];
}
arr1v(rho0,i) = arr1v(rho0,i) + rhoa0j;
arr1v(rho0,j) = arr1v(rho0,j) + rhoa0i;
// For ialloy = 2, use single-element value (not average)
if (meam_data.ialloy!=2) {
arr2v(t_ave,1,i) = arr2v(t_ave,1,i) + meam_data.t1_meam[eltj]*rhoa0j;
arr2v(t_ave,2,i) = arr2v(t_ave,2,i) + meam_data.t2_meam[eltj]*rhoa0j;
arr2v(t_ave,3,i) = arr2v(t_ave,3,i) + meam_data.t3_meam[eltj]*rhoa0j;
arr2v(t_ave,1,j) = arr2v(t_ave,1,j) + meam_data.t1_meam[elti]*rhoa0i;
arr2v(t_ave,2,j) = arr2v(t_ave,2,j) + meam_data.t2_meam[elti]*rhoa0i;
arr2v(t_ave,3,j) = arr2v(t_ave,3,j) + meam_data.t3_meam[elti]*rhoa0i;
}
if (meam_data.ialloy==1) {
arr2v(tsq_ave,1,i) = arr2v(tsq_ave,1,i) + meam_data.t1_meam[eltj]*meam_data.t1_meam[eltj]*rhoa0j;
arr2v(tsq_ave,2,i) = arr2v(tsq_ave,2,i) + meam_data.t2_meam[eltj]*meam_data.t2_meam[eltj]*rhoa0j;
arr2v(tsq_ave,3,i) = arr2v(tsq_ave,3,i) + meam_data.t3_meam[eltj]*meam_data.t3_meam[eltj]*rhoa0j;
arr2v(tsq_ave,1,j) = arr2v(tsq_ave,1,j) + meam_data.t1_meam[elti]*meam_data.t1_meam[elti]*rhoa0i;
arr2v(tsq_ave,2,j) = arr2v(tsq_ave,2,j) + meam_data.t2_meam[elti]*meam_data.t2_meam[elti]*rhoa0i;
arr2v(tsq_ave,3,j) = arr2v(tsq_ave,3,j) + meam_data.t3_meam[elti]*meam_data.t3_meam[elti]*rhoa0i;
}
arr1v(arho2b,i) = arr1v(arho2b,i) + rhoa2j;
arr1v(arho2b,j) = arr1v(arho2b,j) + rhoa2i;
A1j = rhoa1j/rij;
A2j = rhoa2j/rij2;
A3j = rhoa3j/(rij2*rij);
A1i = rhoa1i/rij;
A2i = rhoa2i/rij2;
A3i = rhoa3i/(rij2*rij);
nv2 = 1;
nv3 = 1;
for(m=1; m<=3; m++) {
arr2v(arho1,m,i) = arr2v(arho1,m,i) + A1j*delij[m];
arr2v(arho1,m,j) = arr2v(arho1,m,j) - A1i*delij[m];
arr2v(arho3b,m,i) = arr2v(arho3b,m,i) + rhoa3j*delij[m]/rij;
arr2v(arho3b,m,j) = arr2v(arho3b,m,j) - rhoa3i*delij[m]/rij;
for(n=m; n<=3; n++) {
arr2v(arho2,nv2,i) = arr2v(arho2,nv2,i) + A2j*delij[m]*delij[n];
arr2v(arho2,nv2,j) = arr2v(arho2,nv2,j) + A2i*delij[m]*delij[n];
nv2 = nv2+1;
for(p=n; p<=3; p++) {
arr2v(arho3,nv3,i) = arr2v(arho3,nv3,i) + A3j*delij[m]*delij[n]*delij[p];
arr2v(arho3,nv3,j) = arr2v(arho3,nv3,j) - A3i*delij[m]*delij[n]*delij[p];
nv3 = nv3+1;
}
}
A1j = rhoa1j / rij;
A2j = rhoa2j / rij2;
A3j = rhoa3j / (rij2 * rij);
A1i = rhoa1i / rij;
A2i = rhoa2i / rij2;
A3i = rhoa3i / (rij2 * rij);
nv2 = 1;
nv3 = 1;
for (m = 1; m <= 3; m++) {
arr2v(arho1, m, i) = arr2v(arho1, m, i) + A1j * delij[m];
arr2v(arho1, m, j) = arr2v(arho1, m, j) - A1i * delij[m];
arr2v(arho3b, m, i) = arr2v(arho3b, m, i) + rhoa3j * delij[m] / rij;
arr2v(arho3b, m, j) = arr2v(arho3b, m, j) - rhoa3i * delij[m] / rij;
for (n = m; n <= 3; n++) {
arr2v(arho2, nv2, i) =
arr2v(arho2, nv2, i) + A2j * delij[m] * delij[n];
arr2v(arho2, nv2, j) =
arr2v(arho2, nv2, j) + A2i * delij[m] * delij[n];
nv2 = nv2 + 1;
for (p = n; p <= 3; p++) {
arr2v(arho3, nv3, i) =
arr2v(arho3, nv3, i) + A3j * delij[m] * delij[n] * delij[p];
arr2v(arho3, nv3, j) =
arr2v(arho3, nv3, j) - A3i * delij[m] * delij[n] * delij[p];
nv3 = nv3 + 1;
}
}
}
}
}
}
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void screen(int i, int j, int nmax, double *x, double rijsq, double *sij, int numneigh_full, int *firstneigh_full, int ntype, int *type, int *fmap)
void
screen(int i, int j, int nmax, double* x, double rijsq, double* sij,
int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap)
// Screening function
// Inputs: i = atom 1 id (integer)
// j = atom 2 id (integer)
// rijsq = squared distance between i and j
// Outputs: sij = screening function
{
arrdim2v(x,3,nmax)
{
int k,nk/*,m*/;
int elti,eltj,eltk;
double delxik,delyik,delzik;
double delxjk,delyjk,delzjk;
double riksq,rjksq,xik,xjk,cikj,a,delc,sikj/*,fcij,rij*/;
double Cmax,Cmin,rbound;
arrdim2v(x, 3, nmax)
*sij = 1.0;
eltj = arr1v(fmap,arr1v(type,j));
elti = arr1v(fmap,arr1v(type,j));
int k,
nk /*,m*/;
int elti, eltj, eltk;
double delxik, delyik, delzik;
double delxjk, delyjk, delzjk;
double riksq, rjksq, xik, xjk, cikj, a, delc, sikj /*,fcij,rij*/;
double Cmax, Cmin, rbound;
// if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
rbound = meam_data.ebound_meam[elti][eltj]*rijsq;
*sij = 1.0;
eltj = arr1v(fmap, arr1v(type, j));
elti = arr1v(fmap, arr1v(type, j));
for(nk=1; nk<=numneigh_full; nk++) {
k = arr1v(firstneigh_full,nk);
eltk = arr1v(fmap,arr1v(type,k));
if (k==j) continue;
delxjk = arr2v(x,1,k) - arr2v(x,1,j);
delyjk = arr2v(x,2,k) - arr2v(x,2,j);
delzjk = arr2v(x,3,k) - arr2v(x,3,j);
rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk;
if (rjksq > rbound) continue;
delxik = arr2v(x,1,k) - arr2v(x,1,i);
delyik = arr2v(x,2,k) - arr2v(x,2,i);
delzik = arr2v(x,3,k) - arr2v(x,3,i);
riksq = delxik*delxik + delyik*delyik + delzik*delzik;
if (riksq > rbound) continue;
xik = riksq/rijsq;
xjk = rjksq/rijsq;
a = 1 - (xik-xjk)*(xik-xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a<=0.0) continue;
cikj = (2.0*(xik+xjk) + a - 2.0)/a;
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
if (cikj>=Cmax)
continue;
// note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case here
// (other negative cikj cases were handled by the test on "a" above)
else if (cikj<=Cmin) {
*sij = 0.0;
break;
} else {
delc = Cmax - Cmin;
cikj = (cikj-Cmin)/delc;
fcut(cikj,&sikj);
}
*sij = *sij * sikj;
}
// if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
rbound = meam_data.ebound_meam[elti][eltj] * rijsq;
}
for (nk = 1; nk <= numneigh_full; nk++) {
k = arr1v(firstneigh_full, nk);
eltk = arr1v(fmap, arr1v(type, k));
if (k == j)
continue;
delxjk = arr2v(x, 1, k) - arr2v(x, 1, j);
delyjk = arr2v(x, 2, k) - arr2v(x, 2, j);
delzjk = arr2v(x, 3, k) - arr2v(x, 3, j);
rjksq = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
if (rjksq > rbound)
continue;
delxik = arr2v(x, 1, k) - arr2v(x, 1, i);
delyik = arr2v(x, 2, k) - arr2v(x, 2, i);
delzik = arr2v(x, 3, k) - arr2v(x, 3, i);
riksq = delxik * delxik + delyik * delyik + delzik * delzik;
if (riksq > rbound)
continue;
xik = riksq / rijsq;
xjk = rjksq / rijsq;
a = 1 - (xik - xjk) * (xik - xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a <= 0.0)
continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax)
continue;
// note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case here
// (other negative cikj cases were handled by the test on "a" above)
else if (cikj <= Cmin) {
*sij = 0.0;
break;
} else {
delc = Cmax - Cmin;
cikj = (cikj - Cmin) / delc;
fcut(cikj, &sikj);
}
*sij = *sij * sikj;
}
}
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair)
{
// Inputs: i,j,k = id's of 3 atom triplet
// jn = id of i-j pair
// rij2 = squared distance between i and j
// Outputs: dsij1 = deriv. of sij w.r.t. rik
// dsij2 = deriv. of sij w.r.t. rjk
void
dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2,
double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double* x,
double* scrfcn, double* fcpair)
{
// Inputs: i,j,k = id's of 3 atom triplet
// jn = id of i-j pair
// rij2 = squared distance between i and j
// Outputs: dsij1 = deriv. of sij w.r.t. rik
// dsij2 = deriv. of sij w.r.t. rjk
arrdim2v(x,3,nmax)
int elti,eltj,eltk;
double rik2,rjk2;
arrdim2v(x, 3, nmax) int elti, eltj, eltk;
double rik2, rjk2;
double dxik,dyik,dzik;
double dxjk,dyjk,dzjk;
double rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a;
double Cmax,Cmin,dCikj1,dCikj2;
double dxik, dyik, dzik;
double dxjk, dyjk, dzjk;
double rbound, delc, sij, xik, xjk, cikj, sikj, dfc, a;
double Cmax, Cmin, dCikj1, dCikj2;
sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn);
elti = arr1v(fmap,arr1v(type,i));
eltj = arr1v(fmap,arr1v(type,j));
eltk = arr1v(fmap,arr1v(type,k));
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn);
elti = arr1v(fmap, arr1v(type, i));
eltj = arr1v(fmap, arr1v(type, j));
eltk = arr1v(fmap, arr1v(type, k));
Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
*dsij1 = 0.0;
*dsij2 = 0.0;
if (!iszero(sij) && !iszero(sij-1.0)) {
rbound = rij2*meam_data.ebound_meam[elti][eltj];
delc = Cmax-Cmin;
dxjk = arr2v(x,1,k) - arr2v(x,1,j);
dyjk = arr2v(x,2,k) - arr2v(x,2,j);
dzjk = arr2v(x,3,k) - arr2v(x,3,j);
rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk;
if (rjk2<=rbound) {
dxik = arr2v(x,1,k) - arr2v(x,1,i);
dyik = arr2v(x,2,k) - arr2v(x,2,i);
dzik = arr2v(x,3,k) - arr2v(x,3,i);
rik2 = dxik*dxik + dyik*dyik + dzik*dzik;
if (rik2<=rbound) {
xik = rik2/rij2;
xjk = rjk2/rij2;
a = 1 - (xik-xjk)*(xik-xjk);
if (!iszero(a)) {
cikj = (2.0*(xik+xjk) + a - 2.0)/a;
if (cikj>=Cmin && cikj<=Cmax) {
cikj = (cikj-Cmin)/delc;
dfcut(cikj,&sikj,&dfc);
dCfunc2(rij2,rik2,rjk2,&dCikj1,&dCikj2);
a = sij/delc*dfc/sikj;
*dsij1 = a*dCikj1;
*dsij2 = a*dCikj2;
}
}
*dsij1 = 0.0;
*dsij2 = 0.0;
if (!iszero(sij) && !iszero(sij - 1.0)) {
rbound = rij2 * meam_data.ebound_meam[elti][eltj];
delc = Cmax - Cmin;
dxjk = arr2v(x, 1, k) - arr2v(x, 1, j);
dyjk = arr2v(x, 2, k) - arr2v(x, 2, j);
dzjk = arr2v(x, 3, k) - arr2v(x, 3, j);
rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
if (rjk2 <= rbound) {
dxik = arr2v(x, 1, k) - arr2v(x, 1, i);
dyik = arr2v(x, 2, k) - arr2v(x, 2, i);
dzik = arr2v(x, 3, k) - arr2v(x, 3, i);
rik2 = dxik * dxik + dyik * dyik + dzik * dzik;
if (rik2 <= rbound) {
xik = rik2 / rij2;
xjk = rjk2 / rij2;
a = 1 - (xik - xjk) * (xik - xjk);
if (!iszero(a)) {
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
if (cikj >= Cmin && cikj <= Cmax) {
cikj = (cikj - Cmin) / delc;
dfcut(cikj, &sikj, &dfc);
dCfunc2(rij2, rik2, rjk2, &dCikj1, &dCikj2);
a = sij / delc * dfc / sikj;
*dsij1 = a * dCikj1;
*dsij2 = a * dCikj2;
}
}
}
}
}
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
fcut(double xi, double* fc)
{
// cutoff function
double a;
if (xi >= 1.0)
*fc = 1.0;
else if (xi <= 0.0)
*fc = 0.0;
else {
a = 1.0 - xi;
a = a * a;
a = a * a;
a = 1.0 - a;
*fc = a * a;
// fc = xi
}
}
void fcut(double xi, double * fc)
{
// cutoff function
double a;
if (xi>=1.0)
*fc = 1.0;
else if (xi<=0.0)
*fc = 0.0;
else {
a = 1.0-xi;
a = a*a;
a = a*a;
a = 1.0-a;
*fc = a*a;
// fc = xi
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
dfcut(double xi, double* fc, double* dfc)
{
// cutoff function and its derivative
double a, a3, a4;
if (xi >= 1.0) {
*fc = 1.0;
*dfc = 0.0;
} else if (xi <= 0.0) {
*fc = 0.0;
*dfc = 0.0;
} else {
a = 1.0 - xi;
a3 = a * a * a;
a4 = a * a3;
*fc = pow((1.0 - a4), 2);
*dfc = 8 * (1.0 - a4) * a3;
// fc = xi
// dfc = 1.d0
}
}
void dfcut(double xi, double *fc, double *dfc)
{
// cutoff function and its derivative
double a,a3,a4;
if (xi>=1.0) {
*fc = 1.0;
*dfc = 0.0;
} else if (xi<=0.0) {
*fc = 0.0;
*dfc = 0.0;
} else {
a = 1.0-xi;
a3 = a*a*a;
a4 = a*a3;
*fc = pow((1.0-a4), 2);
*dfc = 8*(1.0-a4)*a3;
// fc = xi
// dfc = 1.d0
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
dCfunc(double rij2, double rik2, double rjk2, double* dCikj)
{
// Inputs: rij,rij2,rik2,rjk2
// Outputs: dCikj = derivative of Cikj w.r.t. rij
double rij4, a, b, denom;
void dCfunc(double rij2,double rik2,double rjk2,double *dCikj)
{
// Inputs: rij,rij2,rik2,rjk2
// Outputs: dCikj = derivative of Cikj w.r.t. rij
double rij4,a,b,denom;
rij4 = rij2 * rij2;
a = rik2 - rjk2;
b = rik2 + rjk2;
denom = rij4 - a * a;
denom = denom * denom;
*dCikj = -4 * (-2 * rij2 * a * a + rij4 * b + a * a * b) / denom;
}
rij4 = rij2*rij2;
a = rik2-rjk2;
b = rik2+rjk2;
denom = rij4 - a*a;
denom = denom*denom;
*dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom;
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCikj2)
{
// Inputs: rij,rij2,rik2,rjk2
// Outputs: dCikj1 = derivative of Cikj w.r.t. rik
// dCikj2 = derivative of Cikj w.r.t. rjk
double rij4, rik4, rjk4, a, b, denom;
void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2)
{
// Inputs: rij,rij2,rik2,rjk2
// Outputs: dCikj1 = derivative of Cikj w.r.t. rik
// dCikj2 = derivative of Cikj w.r.t. rjk
double rij4,rik4,rjk4,a,b,denom;
rij4 = rij2*rij2;
rik4 = rik2*rik2;
rjk4 = rjk2*rjk2;
a = rik2-rjk2;
b = rik2+rjk2;
denom = rij4 - a*a;
denom = denom*denom;
*dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/ denom;
*dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/ denom;
(void)(b);
}
//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
rij4 = rij2 * rij2;
rik4 = rik2 * rik2;
rjk4 = rjk2 * rjk2;
a = rik2 - rjk2;
b = rik2 + rjk2;
denom = rij4 - a * a;
denom = denom * denom;
*dCikj1 = 4 * rij2 *
(rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;
*dCikj2 = 4 * rij2 *
(rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;
(void)(b);
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
}

1077
src/USER-MEAMC/meam_force.cpp
View File

File diff suppressed because it is too large Load Diff

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

File diff suppressed because it is too large Load Diff

148
src/USER-MEAMC/meam_setup_global.cpp
View File

@ -17,85 +17,83 @@ extern "C" {
//
//
void meam_setup_global_(int *nelt, int *lat, double *z, int *ielement, double *atwt, double *alpha,
double *b0, double *b1, double *b2, double *b3, double *alat, double *esub, double *asub,
double *t0, double *t1, double *t2, double *t3, double *rozero, int *ibar)
{
int i;
double tmplat[maxelt];
void
meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt,
double* alpha, double* b0, double* b1, double* b2,
double* b3, double* alat, double* esub, double* asub,
double* t0, double* t1, double* t2, double* t3,
double* rozero, int* ibar)
{
meam_data.neltypes = *nelt;
int i;
double tmplat[maxelt];
for(i=1; i<=*nelt; i++) {
if (arr1v(lat,i)==0)
meam_data.lattce_meam[i][i] = FCC;
else if (arr1v(lat,i)==1)
meam_data.lattce_meam[i][i] = BCC;
else if (arr1v(lat,i)==2)
meam_data.lattce_meam[i][i] = HCP;
else if (arr1v(lat,i)==3)
meam_data.lattce_meam[i][i] = DIM;
else if (arr1v(lat,i)==4)
meam_data.lattce_meam[i][i] = DIA;
else {
// unknown
}
meam_data.neltypes = *nelt;
meam_data.Z_meam[i] = arr1v(z,i);
meam_data.ielt_meam[i] = arr1v(ielement,i);
meam_data.alpha_meam[i][i] = arr1v(alpha,i);
meam_data.beta0_meam[i] = arr1v(b0,i);
meam_data.beta1_meam[i] = arr1v(b1,i);
meam_data.beta2_meam[i] = arr1v(b2,i);
meam_data.beta3_meam[i] = arr1v(b3,i);
tmplat[i] = arr1v(alat,i);
meam_data.Ec_meam[i][i] = arr1v(esub,i);
meam_data.A_meam[i] = arr1v(asub,i);
meam_data.t0_meam[i] = arr1v(t0,i);
meam_data.t1_meam[i] = arr1v(t1,i);
meam_data.t2_meam[i] = arr1v(t2,i);
meam_data.t3_meam[i] = arr1v(t3,i);
meam_data.rho0_meam[i] = arr1v(rozero,i);
meam_data.ibar_meam[i] = arr1v(ibar,i);
for (i = 1; i <= *nelt; i++) {
if (arr1v(lat, i) == 0)
meam_data.lattce_meam[i][i] = FCC;
else if (arr1v(lat, i) == 1)
meam_data.lattce_meam[i][i] = BCC;
else if (arr1v(lat, i) == 2)
meam_data.lattce_meam[i][i] = HCP;
else if (arr1v(lat, i) == 3)
meam_data.lattce_meam[i][i] = DIM;
else if (arr1v(lat, i) == 4)
meam_data.lattce_meam[i][i] = DIA;
else {
// unknown
}
if (meam_data.lattce_meam[i][i]==FCC)
meam_data.re_meam[i][i] = tmplat[i]/sqrt(2.0);
else if (meam_data.lattce_meam[i][i]==BCC)
meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/2.0;
else if (meam_data.lattce_meam[i][i]==HCP)
meam_data.re_meam[i][i] = tmplat[i];
else if (meam_data.lattce_meam[i][i]==DIM)
meam_data.re_meam[i][i] = tmplat[i];
else if (meam_data.lattce_meam[i][i]==DIA)
meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/4.0;
else {
// error
}
}
// Set some defaults
meam_data.rc_meam = 4.0;
meam_data.delr_meam = 0.1;
setall2d(meam_data.attrac_meam, 0.0);
setall2d(meam_data.repuls_meam, 0.0);
setall3d(meam_data.Cmax_meam, 2.8);
setall3d(meam_data.Cmin_meam, 2.0);
setall2d(meam_data.ebound_meam, pow(2.8,2)/(4.0*(2.8-1.0)));
setall2d(meam_data.delta_meam, 0.0);
setall2d(meam_data.nn2_meam, 0);
setall2d(meam_data.zbl_meam, 1);
meam_data.gsmooth_factor = 99.0;
meam_data.augt1 = 1;
meam_data.ialloy = 0;
meam_data.mix_ref_t = 0;
meam_data.emb_lin_neg = 0;
meam_data.bkgd_dyn = 0;
meam_data.erose_form = 0;
}
meam_data.Z_meam[i] = arr1v(z, i);
meam_data.ielt_meam[i] = arr1v(ielement, i);
meam_data.alpha_meam[i][i] = arr1v(alpha, i);
meam_data.beta0_meam[i] = arr1v(b0, i);
meam_data.beta1_meam[i] = arr1v(b1, i);
meam_data.beta2_meam[i] = arr1v(b2, i);
meam_data.beta3_meam[i] = arr1v(b3, i);
tmplat[i] = arr1v(alat, i);
meam_data.Ec_meam[i][i] = arr1v(esub, i);
meam_data.A_meam[i] = arr1v(asub, i);
meam_data.t0_meam[i] = arr1v(t0, i);
meam_data.t1_meam[i] = arr1v(t1, i);
meam_data.t2_meam[i] = arr1v(t2, i);
meam_data.t3_meam[i] = arr1v(t3, i);
meam_data.rho0_meam[i] = arr1v(rozero, i);
meam_data.ibar_meam[i] = arr1v(ibar, i);
if (meam_data.lattce_meam[i][i] == FCC)
meam_data.re_meam[i][i] = tmplat[i] / sqrt(2.0);
else if (meam_data.lattce_meam[i][i] == BCC)
meam_data.re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
else if (meam_data.lattce_meam[i][i] == HCP)
meam_data.re_meam[i][i] = tmplat[i];
else if (meam_data.lattce_meam[i][i] == DIM)
meam_data.re_meam[i][i] = tmplat[i];
else if (meam_data.lattce_meam[i][i] == DIA)
meam_data.re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
else {
// error
}
}
// Set some defaults
meam_data.rc_meam = 4.0;
meam_data.delr_meam = 0.1;
setall2d(meam_data.attrac_meam, 0.0);
setall2d(meam_data.repuls_meam, 0.0);
setall3d(meam_data.Cmax_meam, 2.8);
setall3d(meam_data.Cmin_meam, 2.0);
setall2d(meam_data.ebound_meam, pow(2.8, 2) / (4.0 * (2.8 - 1.0)));
setall2d(meam_data.delta_meam, 0.0);
setall2d(meam_data.nn2_meam, 0);
setall2d(meam_data.zbl_meam, 1);
meam_data.gsmooth_factor = 99.0;
meam_data.augt1 = 1;
meam_data.ialloy = 0;
meam_data.mix_ref_t = 0;
meam_data.emb_lin_neg = 0;
meam_data.bkgd_dyn = 0;
meam_data.erose_form = 0;
}
}

332
src/USER-MEAMC/meam_setup_param.cpp
View File

@ -3,22 +3,23 @@ extern "C" {
//
// do a sanity check on index parameters
void meam_checkindex(int num, int lim, int nidx, int *idx /*idx(3)*/, int *ierr)
{
//: idx[0..2]
*ierr = 0;
if (nidx < num) {
*ierr = 2;
return;
}
void
meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr)
{
//: idx[0..2]
*ierr = 0;
if (nidx < num) {
*ierr = 2;
return;
}
for (int i=0; i<num; i++) {
if ((idx[i] < 1) || (idx[i] > lim)) {
*ierr = 3;
return;
}
}
}
for (int i = 0; i < num; i++) {
if ((idx[i] < 1) || (idx[i] > lim)) {
*ierr = 3;
return;
}
}
}
//
// Declaration in pair_meam.h:
@ -56,171 +57,184 @@ extern "C" {
// 19 = emb_lin_neg
// 20 = bkgd_dyn
void meam_setup_param_(int *which_p, double *value_p, int *nindex_p, int *index /*index(3)*/, int *errorflag)
{
//: index[0..2]
int i1, i2, val;
*errorflag = 0;
int which = *which_p;
double value = *value_p;
int nindex = *nindex_p;
void
meam_setup_param_(int* which_p, double* value_p, int* nindex_p,
int* index /*index(3)*/, int* errorflag)
{
//: index[0..2]
int i1, i2, val;
*errorflag = 0;
int which = *which_p;
double value = *value_p;
int nindex = *nindex_p;
switch(which) {
// 0 = Ec_meam
case 0:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.Ec_meam[index[0]][index[1]] = value;
break;
switch (which) {
// 0 = Ec_meam
case 0:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.Ec_meam[index[0]][index[1]] = value;
break;
// 1 = alpha_meam
case 1:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.alpha_meam[index[0]][index[1]] = value;
break;
// 1 = alpha_meam
case 1:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.alpha_meam[index[0]][index[1]] = value;
break;
// 2 = rho0_meam
case 2:
meam_checkindex(1,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.rho0_meam[index[0]] = value;
break;
// 2 = rho0_meam
case 2:
meam_checkindex(1, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.rho0_meam[index[0]] = value;
break;
// 3 = delta_meam
case 3:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.delta_meam[index[0]][index[1]] = value;
break;
// 3 = delta_meam
case 3:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.delta_meam[index[0]][index[1]] = value;
break;
// 4 = lattce_meam
case 4:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
val = (int)value;
// 4 = lattce_meam
case 4:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
val = (int)value;
if (val==0)
meam_data.lattce_meam[index[0]][index[1]] = FCC;
else if (val==1)
meam_data.lattce_meam[index[0]][index[1]] = BCC;
else if (val==2)
meam_data.lattce_meam[index[0]][index[1]] = HCP;
else if (val==3)
meam_data.lattce_meam[index[0]][index[1]] = DIM;
else if (val==4)
meam_data.lattce_meam[index[0]][index[1]] = DIA;
else if (val==5)
meam_data.lattce_meam[index[0]][index[1]] = B1;
else if (val==6)
meam_data.lattce_meam[index[0]][index[1]] = C11;
else if (val==7)
meam_data.lattce_meam[index[0]][index[1]] = L12;
else if (val==8)
meam_data.lattce_meam[index[0]][index[1]] = B2;
break;
if (val == 0)
meam_data.lattce_meam[index[0]][index[1]] = FCC;
else if (val == 1)
meam_data.lattce_meam[index[0]][index[1]] = BCC;
else if (val == 2)
meam_data.lattce_meam[index[0]][index[1]] = HCP;
else if (val == 3)
meam_data.lattce_meam[index[0]][index[1]] = DIM;
else if (val == 4)
meam_data.lattce_meam[index[0]][index[1]] = DIA;
else if (val == 5)
meam_data.lattce_meam[index[0]][index[1]] = B1;
else if (val == 6)
meam_data.lattce_meam[index[0]][index[1]] = C11;
else if (val == 7)
meam_data.lattce_meam[index[0]][index[1]] = L12;
else if (val == 8)
meam_data.lattce_meam[index[0]][index[1]] = B2;
break;
// 5 = attrac_meam
case 5:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.attrac_meam[index[0]][index[1]] = value;
break;
// 5 = attrac_meam
case 5:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.attrac_meam[index[0]][index[1]] = value;
break;
// 6 = repuls_meam
case 6:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.repuls_meam[index[0]][index[1]] = value;
break;
// 6 = repuls_meam
case 6:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.repuls_meam[index[0]][index[1]] = value;
break;
// 7 = nn2_meam
case 7:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
i1 = min(index[0],index[1]);
i2 = max(index[0],index[1]);
meam_data.nn2_meam[i1][i2] = (int)value;
break;
// 7 = nn2_meam
case 7:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
i1 = min(index[0], index[1]);
i2 = max(index[0], index[1]);
meam_data.nn2_meam[i1][i2] = (int)value;
break;
// 8 = Cmin_meam
case 8:
meam_checkindex(3,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value;
break;
// 8 = Cmin_meam
case 8:
meam_checkindex(3, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value;
break;
// 9 = Cmax_meam
case 9:
meam_checkindex(3,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value;
break;
// 9 = Cmax_meam
case 9:
meam_checkindex(3, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value;
break;
// 10 = rc_meam
case 10:
meam_data.rc_meam = value;
break;
// 10 = rc_meam
case 10:
meam_data.rc_meam = value;
break;
// 11 = delr_meam
case 11:
meam_data.delr_meam = value;
break;
// 11 = delr_meam
case 11:
meam_data.delr_meam = value;
break;
// 12 = augt1
case 12:
meam_data.augt1 = (int)value;
break;
// 12 = augt1
case 12:
meam_data.augt1 = (int)value;
break;
// 13 = gsmooth
case 13:
meam_data.gsmooth_factor = value;
break;
// 13 = gsmooth
case 13:
meam_data.gsmooth_factor = value;
break;
// 14 = re_meam
case 14:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
meam_data.re_meam[index[0]][index[1]] = value;
break;
// 14 = re_meam
case 14:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
meam_data.re_meam[index[0]][index[1]] = value;
break;
// 15 = ialloy
case 15:
meam_data.ialloy = (int)value;
break;
// 15 = ialloy
case 15:
meam_data.ialloy = (int)value;
break;
// 16 = mixture_ref_t
case 16:
meam_data.mix_ref_t = (int)value;
break;
// 16 = mixture_ref_t
case 16:
meam_data.mix_ref_t = (int)value;
break;
// 17 = erose_form
case 17:
meam_data.erose_form = (int)value;
break;
// 17 = erose_form
case 17:
meam_data.erose_form = (int)value;
break;
// 18 = zbl_meam
case 18:
meam_checkindex(2,maxelt,nindex,index,errorflag);
if (*errorflag!=0) return;
i1 = min(index[0],index[1]);
i2 = max(index[0],index[1]);
meam_data.zbl_meam[i1][i2] = (int)value;
break;
// 18 = zbl_meam
case 18:
meam_checkindex(2, maxelt, nindex, index, errorflag);
if (*errorflag != 0)
return;
i1 = min(index[0], index[1]);
i2 = max(index[0], index[1]);
meam_data.zbl_meam[i1][i2] = (int)value;
break;
// 19 = emb_lin_neg
case 19:
meam_data.emb_lin_neg = (int)value;
break;
// 19 = emb_lin_neg
case 19:
meam_data.emb_lin_neg = (int)value;
break;
// 20 = bkgd_dyn
case 20:
meam_data.bkgd_dyn = (int)value;
break;
default:
*errorflag = 1;
}
}
// 20 = bkgd_dyn
case 20:
meam_data.bkgd_dyn = (int)value;
break;
default:
*errorflag = 1;
}
}
}