255 lines
9.1 KiB
C
255 lines
9.1 KiB
C
|
#include <stdio.h>
|
||
|
|
#include <stdlib.h>
|
||
|
|
|
||
|
|
#include <allocate.h>
|
||
|
|
#include <atom.h>
|
||
|
|
#include <eam.h>
|
||
|
|
#include <parameter.h>
|
||
|
|
|
||
|
|
void initEam(Eam* eam, const char *input_file, int ntypes) {
|
||
|
|
eam->nmax = 0;
|
||
|
|
eam->fp = NULL;
|
||
|
|
eam->ntypes = ntypes;
|
||
|
|
eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT));
|
||
|
|
}
|
||
|
|
|
||
|
|
void coeff(Eam* eam, const char* arg) {
|
||
|
|
read_file(eam->file, arg);
|
||
|
|
int n = strlen(arg) + 1;
|
||
|
|
int ntypes = eam->ntypes;
|
||
|
|
double cutmax = eam->file->cut;
|
||
|
|
for(int i=0; i<ntypes*ntypes; i++)
|
||
|
|
eam->cutforcesq[i] = cutmax * cutmax;
|
||
|
|
}
|
||
|
|
|
||
|
|
void init_style(Eam* eam) {
|
||
|
|
// convert read-in file(s) to arrays and spline them
|
||
|
|
file2array(eam);
|
||
|
|
array2spline(eam);
|
||
|
|
}
|
||
|
|
|
||
|
|
void read_file(Funcfl* file, const char* filename) {
|
||
|
|
FILE* fptr;
|
||
|
|
char line[MAXLINE];
|
||
|
|
|
||
|
|
fptr = fopen(filename, "r");
|
||
|
|
if(fptr == NULL) {
|
||
|
|
printf("Can't open EAM Potential file: %s\n", filename);
|
||
|
|
exit(0);
|
||
|
|
}
|
||
|
|
|
||
|
|
int tmp;
|
||
|
|
fgets(line, MAXLINE, fptr);
|
||
|
|
fgets(line, MAXLINE, fptr);
|
||
|
|
sscanf(line, "%d %lg", &tmp, &(file->mass));
|
||
|
|
fgets(line, MAXLINE, fptr);
|
||
|
|
sscanf(line, "%d %lg %d %lg %lg", &file->nrho, &file->drho, &file->nr, &file->dr, &file->cut);
|
||
|
|
|
||
|
|
//printf("Read: %lf %i %lf %i %lf %lf\n",file->mass,file->nrho,file->drho,file->nr,file->dr,file->cut);
|
||
|
|
file->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nhro + 1) * sizeof(MD_FLOAT));
|
||
|
|
file->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT));
|
||
|
|
file->zr = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT));
|
||
|
|
grab(fptr, file->nrho, file->frho);
|
||
|
|
grab(fptr, file->nr, file->zr);
|
||
|
|
grab(fptr, file->nr, file->rhor);
|
||
|
|
for(int i = file->nrho; i > 0; i--) file->frho[i] = file->frho[i - 1];
|
||
|
|
for(int i = file->nr; i > 0; i--) file->rhor[i] = file->rhor[i - 1];
|
||
|
|
for(int i = file->nr; i > 0; i--) file->zr[i] = file->zr[i - 1];
|
||
|
|
fclose(fptr);
|
||
|
|
}
|
||
|
|
|
||
|
|
void file2array(Eam* eam) {
|
||
|
|
int i, j, k, m, n;
|
||
|
|
double sixth = 1.0 / 6.0;
|
||
|
|
|
||
|
|
// determine max function params from all active funcfl files
|
||
|
|
// active means some element is pointing at it via map
|
||
|
|
int active;
|
||
|
|
double rmax, rhomax;
|
||
|
|
dr = drho = rmax = rhomax = 0.0;
|
||
|
|
active = 0;
|
||
|
|
Funcfl* file = eam->file;
|
||
|
|
dr = MAX(dr, file->dr);
|
||
|
|
drho = MAX(drho, file->drho);
|
||
|
|
rmax = MAX(rmax, (file->nr - 1) * file->dr);
|
||
|
|
rhomax = MAX(rhomax, (file->nrho - 1) * file->drho);
|
||
|
|
|
||
|
|
// set nr,nrho from cutoff and spacings
|
||
|
|
// 0.5 is for round-off in divide
|
||
|
|
eam->nr = static_cast<int>(rmax / dr + 0.5);
|
||
|
|
eam->nrho = static_cast<int>(rhomax / drho + 0.5);
|
||
|
|
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
// setup frho arrays
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
|
||
|
|
// allocate frho arrays
|
||
|
|
// nfrho = # of funcfl files + 1 for zero array
|
||
|
|
eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nhro + 1) * sizeof(MD_FLOAT));
|
||
|
|
|
||
|
|
// interpolate each file's frho to a single grid and cutoff
|
||
|
|
double r, p, cof1, cof2, cof3, cof4;
|
||
|
|
n = 0;
|
||
|
|
|
||
|
|
for(m = 1; m <= eam->nrho; m++) {
|
||
|
|
r = (m - 1) * drho;
|
||
|
|
p = r / file->drho + 1.0;
|
||
|
|
k = static_cast<int>(p);
|
||
|
|
k = MIN(k, file->nrho - 2);
|
||
|
|
k = MAX(k, 2);
|
||
|
|
p -= k;
|
||
|
|
p = MIN(p, 2.0);
|
||
|
|
cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
|
||
|
|
cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
|
||
|
|
cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
|
||
|
|
cof4 = sixth * p * (p * p - 1.0);
|
||
|
|
eam->frho[m] = cof1 * file->frho[k - 1] + cof2 * file->frho[k] +
|
||
|
|
cof3 * file->frho[k + 1] + cof4 * file->frho[k + 2];
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
// setup rhor arrays
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
|
||
|
|
// allocate rhor arrays
|
||
|
|
// nrhor = # of funcfl files
|
||
|
|
eam->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT));
|
||
|
|
|
||
|
|
// interpolate each file's rhor to a single grid and cutoff
|
||
|
|
for(m = 1; m <= eam->nr; m++) {
|
||
|
|
r = (m - 1) * dr;
|
||
|
|
p = r / file->dr + 1.0;
|
||
|
|
k = static_cast<int>(p);
|
||
|
|
k = MIN(k, file->nr - 2);
|
||
|
|
k = MAX(k, 2);
|
||
|
|
p -= k;
|
||
|
|
p = MIN(p, 2.0);
|
||
|
|
cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
|
||
|
|
cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
|
||
|
|
cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
|
||
|
|
cof4 = sixth * p * (p * p - 1.0);
|
||
|
|
eam->rhor[m] = cof1 * file->rhor[k - 1] + cof2 * file->rhor[k] +
|
||
|
|
cof3 * file->rhor[k + 1] + cof4 * file->rhor[k + 2];
|
||
|
|
//if(m==119)printf("BuildRho: %e %e %e %e %e %e\n",rhor[m],cof1,cof2,cof3,cof4,file->rhor[k]);
|
||
|
|
}
|
||
|
|
|
||
|
|
// type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to
|
||
|
|
// for funcfl files, I,J mapping only depends on I
|
||
|
|
// OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used
|
||
|
|
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
// setup z2r arrays
|
||
|
|
// ------------------------------------------------------------------
|
||
|
|
|
||
|
|
// allocate z2r arrays
|
||
|
|
// nz2r = N*(N+1)/2 where N = # of funcfl files
|
||
|
|
eam->z2r = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT));
|
||
|
|
|
||
|
|
// create a z2r array for each file against other files, only for I >= J
|
||
|
|
// interpolate zri and zrj to a single grid and cutoff
|
||
|
|
double zri, zrj;
|
||
|
|
Funcfl* ifile = eam->file;
|
||
|
|
Funcfl* jfile = eam->file;
|
||
|
|
|
||
|
|
for(m = 1; m <= eam->nr; m++) {
|
||
|
|
r = (m - 1) * dr;
|
||
|
|
p = r / ifile->dr + 1.0;
|
||
|
|
k = static_cast<int>(p);
|
||
|
|
k = MIN(k, ifile->nr - 2);
|
||
|
|
k = MAX(k, 2);
|
||
|
|
p -= k;
|
||
|
|
p = MIN(p, 2.0);
|
||
|
|
cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
|
||
|
|
cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
|
||
|
|
cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
|
||
|
|
cof4 = sixth * p * (p * p - 1.0);
|
||
|
|
zri = cof1 * ifile->zr[k - 1] + cof2 * ifile->zr[k] +
|
||
|
|
cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2];
|
||
|
|
|
||
|
|
p = r / jfile->dr + 1.0;
|
||
|
|
k = static_cast<int>(p);
|
||
|
|
k = MIN(k, jfile->nr - 2);
|
||
|
|
k = MAX(k, 2);
|
||
|
|
p -= k;
|
||
|
|
p = MIN(p, 2.0);
|
||
|
|
cof1 = -sixth * p * (p - 1.0) * (p - 2.0);
|
||
|
|
cof2 = 0.5 * (p * p - 1.0) * (p - 2.0);
|
||
|
|
cof3 = -0.5 * p * (p + 1.0) * (p - 2.0);
|
||
|
|
cof4 = sixth * p * (p * p - 1.0);
|
||
|
|
zrj = cof1 * jfile->zr[k - 1] + cof2 * jfile->zr[k] +
|
||
|
|
cof3 * jfile->zr[k + 1] + cof4 * jfile->zr[k + 2];
|
||
|
|
|
||
|
|
eam->z2r[m] = 27.2 * 0.529 * zri * zrj;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
void array2spline(Eam* eam) {
|
||
|
|
rdr = 1.0 / eam->dr;
|
||
|
|
rdrho = 1.0 / eam->drho;
|
||
|
|
nrho_tot = (eam->nrho + 1) * 7 + 64;
|
||
|
|
nr_tot = (eam->nr + 1) * 7 + 64;
|
||
|
|
nrho_tot -= nrho_tot%64;
|
||
|
|
nr_tot -= nr_tot%64;
|
||
|
|
|
||
|
|
int ntypes = eam->ntypes;
|
||
|
|
eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nrho_tot * sizeof(MD_FLOAT));
|
||
|
|
eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT));
|
||
|
|
eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT));
|
||
|
|
interpolate(eam->nrho, eam->drho, eam->frho, eam->frho_spline);
|
||
|
|
interpolate(eam->nr, eam->dr, eam->rhor, eam->rhor_spline);
|
||
|
|
interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline);
|
||
|
|
|
||
|
|
// replicate data for multiple types;
|
||
|
|
for(int tt = 0; tt < ntypes * ntypes; tt++) {
|
||
|
|
for(int k = 0; k < nrho_tot; k++)
|
||
|
|
eam->frho_spline[tt*nrho_tot + k] = eam->frho_spline[k];
|
||
|
|
for(int k = 0; k < nr_tot; k++)
|
||
|
|
eam->rhor_spline[tt*nr_tot + k] = eam->rhor_spline[k];
|
||
|
|
for(int k = 0; k < nr_tot; k++)
|
||
|
|
eam->z2r_spline[tt*nr_tot + k] = eam->z2r_spline[k];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) {
|
||
|
|
for(int m = 1; m <= n; m++) spline[m * 7 + 6] = f[m];
|
||
|
|
|
||
|
|
spline[1 * 7 + 5] = spline[2 * 7 + 6] - spline[1 * 7 + 6];
|
||
|
|
spline[2 * 7 + 5] = 0.5 * (spline[3 * 7 + 6] - spline[1 * 7 + 6]);
|
||
|
|
spline[(n - 1) * 7 + 5] = 0.5 * (spline[n * 7 + 6] - spline[(n - 2) * 7 + 6]);
|
||
|
|
spline[n * 7 + 5] = spline[n * 7 + 6] - spline[(n - 1) * 7 + 6];
|
||
|
|
|
||
|
|
for(int m = 3; m <= n - 2; m++)
|
||
|
|
spline[m * 7 + 5] = ((spline[(m - 2) * 7 + 6] - spline[(m + 2) * 7 + 6]) +
|
||
|
|
8.0 * (spline[(m + 1) * 7 + 6] - spline[(m - 1) * 7 + 6])) / 12.0;
|
||
|
|
|
||
|
|
for(int m = 1; m <= n - 1; m++) {
|
||
|
|
spline[m * 7 + 4] = 3.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]) -
|
||
|
|
2.0 * spline[m * 7 + 5] - spline[(m + 1) * 7 + 5];
|
||
|
|
spline[m * 7 + 3] = spline[m * 7 + 5] + spline[(m + 1) * 7 + 5] -
|
||
|
|
2.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]);
|
||
|
|
}
|
||
|
|
|
||
|
|
spline[n * 7 + 4] = 0.0;
|
||
|
|
spline[n * 7 + 3] = 0.0;
|
||
|
|
|
||
|
|
for(int m = 1; m <= n; m++) {
|
||
|
|
spline[m * 7 + 2] = spline[m * 7 + 5] / delta;
|
||
|
|
spline[m * 7 + 1] = 2.0 * spline[m * 7 + 4] / delta;
|
||
|
|
spline[m * 7 + 0] = 3.0 * spline[m * 7 + 3] / delta;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
void grab(FILE* fptr, MMD_int n, MMD_float* list) {
|
||
|
|
char* ptr;
|
||
|
|
char line[MAXLINE];
|
||
|
|
int i = 0;
|
||
|
|
|
||
|
|
while(i < n) {
|
||
|
|
fgets(line, MAXLINE, fptr);
|
||
|
|
ptr = strtok(line, " \t\n\r\f");
|
||
|
|
list[i++] = atof(ptr);
|
||
|
|
while(ptr = strtok(NULL, " \t\n\r\f")) list[i++] = atof(ptr);
|
||
|
|
}
|
||
|
|
}
|