Fix errors introduced by last changes
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
parent
2dac10469c
commit
40ddc9ad50
@ -41,14 +41,12 @@ void initAtom(Atom *atom)
|
|||||||
atom->Nlocal = 0;
|
atom->Nlocal = 0;
|
||||||
atom->Nghost = 0;
|
atom->Nghost = 0;
|
||||||
atom->Nmax = 0;
|
atom->Nmax = 0;
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->type = NULL;
|
atom->type = NULL;
|
||||||
atom->ntypes = 0;
|
atom->ntypes = 0;
|
||||||
atom->epsilon = NULL;
|
atom->epsilon = NULL;
|
||||||
atom->sigma6 = NULL;
|
atom->sigma6 = NULL;
|
||||||
atom->cutforcesq = NULL;
|
atom->cutforcesq = NULL;
|
||||||
atom->cutneighsq = NULL;
|
atom->cutneighsq = NULL;
|
||||||
#endif
|
|
||||||
}
|
}
|
||||||
|
|
||||||
void createAtom(Atom *atom, Parameter *param)
|
void createAtom(Atom *atom, Parameter *param)
|
||||||
@ -58,8 +56,6 @@ void createAtom(Atom *atom, Parameter *param)
|
|||||||
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd;
|
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd;
|
||||||
atom->Natoms = 4 * param->nx * param->ny * param->nz;
|
atom->Natoms = 4 * param->nx * param->ny * param->nz;
|
||||||
atom->Nlocal = 0;
|
atom->Nlocal = 0;
|
||||||
|
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->ntypes = param->ntypes;
|
atom->ntypes = param->ntypes;
|
||||||
atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
||||||
atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
||||||
@ -71,7 +67,6 @@ void createAtom(Atom *atom, Parameter *param)
|
|||||||
atom->cutneighsq[i] = param->cutneigh * param->cutneigh;
|
atom->cutneighsq[i] = param->cutneigh * param->cutneigh;
|
||||||
atom->cutforcesq[i] = param->cutforce * param->cutforce;
|
atom->cutforcesq[i] = param->cutforce * param->cutforce;
|
||||||
}
|
}
|
||||||
#endif
|
|
||||||
|
|
||||||
MD_FLOAT alat = pow((4.0 / param->rho), (1.0 / 3.0));
|
MD_FLOAT alat = pow((4.0 / param->rho), (1.0 / 3.0));
|
||||||
int ilo = (int) (xlo / (0.5 * alat) - 1);
|
int ilo = (int) (xlo / (0.5 * alat) - 1);
|
||||||
@ -142,9 +137,7 @@ void createAtom(Atom *atom, Parameter *param)
|
|||||||
atom->vx[atom->Nlocal] = vxtmp;
|
atom->vx[atom->Nlocal] = vxtmp;
|
||||||
atom->vy[atom->Nlocal] = vytmp;
|
atom->vy[atom->Nlocal] = vytmp;
|
||||||
atom->vz[atom->Nlocal] = vztmp;
|
atom->vz[atom->Nlocal] = vztmp;
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->type[atom->Nlocal] = rand() % atom->ntypes;
|
atom->type[atom->Nlocal] = rand() % atom->ntypes;
|
||||||
#endif
|
|
||||||
atom->Nlocal++;
|
atom->Nlocal++;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -177,9 +170,7 @@ void growAtom(Atom *atom)
|
|||||||
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||||
atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||||
atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
|
atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
|
||||||
#endif
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
@ -1,16 +1,23 @@
|
|||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
|
||||||
#include <allocate.h>
|
#include <allocate.h>
|
||||||
#include <atom.h>
|
#include <atom.h>
|
||||||
#include <eam.h>
|
#include <eam.h>
|
||||||
#include <parameter.h>
|
#include <parameter.h>
|
||||||
|
#include <util.h>
|
||||||
|
|
||||||
|
#ifndef MAXLINE
|
||||||
|
#define MAXLINE 4096
|
||||||
|
#endif
|
||||||
|
|
||||||
void initEam(Eam* eam, const char *input_file, int ntypes) {
|
void initEam(Eam* eam, const char *input_file, int ntypes) {
|
||||||
eam->nmax = 0;
|
eam->nmax = 0;
|
||||||
eam->fp = NULL;
|
eam->fp = NULL;
|
||||||
eam->ntypes = ntypes;
|
eam->ntypes = ntypes;
|
||||||
eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT));
|
eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT));
|
||||||
|
coeff(eam, input_file);
|
||||||
}
|
}
|
||||||
|
|
||||||
void coeff(Eam* eam, const char* arg) {
|
void coeff(Eam* eam, const char* arg) {
|
||||||
@ -46,7 +53,7 @@ void read_file(Funcfl* file, const char* filename) {
|
|||||||
sscanf(line, "%d %lg %d %lg %lg", &file->nrho, &file->drho, &file->nr, &file->dr, &file->cut);
|
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);
|
//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->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nrho + 1) * sizeof(MD_FLOAT));
|
||||||
file->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 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));
|
file->zr = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT));
|
||||||
grab(fptr, file->nrho, file->frho);
|
grab(fptr, file->nrho, file->frho);
|
||||||
@ -66,18 +73,18 @@ void file2array(Eam* eam) {
|
|||||||
// active means some element is pointing at it via map
|
// active means some element is pointing at it via map
|
||||||
int active;
|
int active;
|
||||||
double rmax, rhomax;
|
double rmax, rhomax;
|
||||||
dr = drho = rmax = rhomax = 0.0;
|
eam->dr = eam->drho = rmax = rhomax = 0.0;
|
||||||
active = 0;
|
active = 0;
|
||||||
Funcfl* file = eam->file;
|
Funcfl* file = eam->file;
|
||||||
dr = MAX(dr, file->dr);
|
eam->dr = MAX(eam->dr, file->dr);
|
||||||
drho = MAX(drho, file->drho);
|
eam->drho = MAX(eam->drho, file->drho);
|
||||||
rmax = MAX(rmax, (file->nr - 1) * file->dr);
|
rmax = MAX(rmax, (file->nr - 1) * file->dr);
|
||||||
rhomax = MAX(rhomax, (file->nrho - 1) * file->drho);
|
rhomax = MAX(rhomax, (file->nrho - 1) * file->drho);
|
||||||
|
|
||||||
// set nr,nrho from cutoff and spacings
|
// set nr,nrho from cutoff and spacings
|
||||||
// 0.5 is for round-off in divide
|
// 0.5 is for round-off in divide
|
||||||
eam->nr = static_cast<int>(rmax / dr + 0.5);
|
eam->nr = (int)(rmax / eam->dr + 0.5);
|
||||||
eam->nrho = static_cast<int>(rhomax / drho + 0.5);
|
eam->nrho = (int)(rhomax / eam->drho + 0.5);
|
||||||
|
|
||||||
// ------------------------------------------------------------------
|
// ------------------------------------------------------------------
|
||||||
// setup frho arrays
|
// setup frho arrays
|
||||||
@ -85,16 +92,16 @@ void file2array(Eam* eam) {
|
|||||||
|
|
||||||
// allocate frho arrays
|
// allocate frho arrays
|
||||||
// nfrho = # of funcfl files + 1 for zero array
|
// nfrho = # of funcfl files + 1 for zero array
|
||||||
eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nhro + 1) * sizeof(MD_FLOAT));
|
eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nrho + 1) * sizeof(MD_FLOAT));
|
||||||
|
|
||||||
// interpolate each file's frho to a single grid and cutoff
|
// interpolate each file's frho to a single grid and cutoff
|
||||||
double r, p, cof1, cof2, cof3, cof4;
|
double r, p, cof1, cof2, cof3, cof4;
|
||||||
n = 0;
|
n = 0;
|
||||||
|
|
||||||
for(m = 1; m <= eam->nrho; m++) {
|
for(m = 1; m <= eam->nrho; m++) {
|
||||||
r = (m - 1) * drho;
|
r = (m - 1) * eam->drho;
|
||||||
p = r / file->drho + 1.0;
|
p = r / file->drho + 1.0;
|
||||||
k = static_cast<int>(p);
|
k = (int)(p);
|
||||||
k = MIN(k, file->nrho - 2);
|
k = MIN(k, file->nrho - 2);
|
||||||
k = MAX(k, 2);
|
k = MAX(k, 2);
|
||||||
p -= k;
|
p -= k;
|
||||||
@ -118,9 +125,9 @@ void file2array(Eam* eam) {
|
|||||||
|
|
||||||
// interpolate each file's rhor to a single grid and cutoff
|
// interpolate each file's rhor to a single grid and cutoff
|
||||||
for(m = 1; m <= eam->nr; m++) {
|
for(m = 1; m <= eam->nr; m++) {
|
||||||
r = (m - 1) * dr;
|
r = (m - 1) * eam->dr;
|
||||||
p = r / file->dr + 1.0;
|
p = r / file->dr + 1.0;
|
||||||
k = static_cast<int>(p);
|
k = (int)(p);
|
||||||
k = MIN(k, file->nr - 2);
|
k = MIN(k, file->nr - 2);
|
||||||
k = MAX(k, 2);
|
k = MAX(k, 2);
|
||||||
p -= k;
|
p -= k;
|
||||||
@ -153,9 +160,9 @@ void file2array(Eam* eam) {
|
|||||||
Funcfl* jfile = eam->file;
|
Funcfl* jfile = eam->file;
|
||||||
|
|
||||||
for(m = 1; m <= eam->nr; m++) {
|
for(m = 1; m <= eam->nr; m++) {
|
||||||
r = (m - 1) * dr;
|
r = (m - 1) * eam->dr;
|
||||||
p = r / ifile->dr + 1.0;
|
p = r / ifile->dr + 1.0;
|
||||||
k = static_cast<int>(p);
|
k = (int)(p);
|
||||||
k = MIN(k, ifile->nr - 2);
|
k = MIN(k, ifile->nr - 2);
|
||||||
k = MAX(k, 2);
|
k = MAX(k, 2);
|
||||||
p -= k;
|
p -= k;
|
||||||
@ -168,7 +175,7 @@ void file2array(Eam* eam) {
|
|||||||
cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2];
|
cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2];
|
||||||
|
|
||||||
p = r / jfile->dr + 1.0;
|
p = r / jfile->dr + 1.0;
|
||||||
k = static_cast<int>(p);
|
k = (int)(p);
|
||||||
k = MIN(k, jfile->nr - 2);
|
k = MIN(k, jfile->nr - 2);
|
||||||
k = MAX(k, 2);
|
k = MAX(k, 2);
|
||||||
p -= k;
|
p -= k;
|
||||||
@ -185,33 +192,33 @@ void file2array(Eam* eam) {
|
|||||||
}
|
}
|
||||||
|
|
||||||
void array2spline(Eam* eam) {
|
void array2spline(Eam* eam) {
|
||||||
rdr = 1.0 / eam->dr;
|
eam->rdr = 1.0 / eam->dr;
|
||||||
rdrho = 1.0 / eam->drho;
|
eam->rdrho = 1.0 / eam->drho;
|
||||||
nrho_tot = (eam->nrho + 1) * 7 + 64;
|
eam->nrho_tot = (eam->nrho + 1) * 7 + 64;
|
||||||
nr_tot = (eam->nr + 1) * 7 + 64;
|
eam->nr_tot = (eam->nr + 1) * 7 + 64;
|
||||||
nrho_tot -= nrho_tot%64;
|
eam->nrho_tot -= eam->nrho_tot%64;
|
||||||
nr_tot -= nr_tot%64;
|
eam->nr_tot -= eam->nr_tot%64;
|
||||||
|
|
||||||
int ntypes = eam->ntypes;
|
int ntypes = eam->ntypes;
|
||||||
eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nrho_tot * sizeof(MD_FLOAT));
|
eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nrho_tot * sizeof(MD_FLOAT));
|
||||||
eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT));
|
eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));
|
||||||
eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT));
|
eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));
|
||||||
interpolate(eam->nrho, eam->drho, eam->frho, eam->frho_spline);
|
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->rhor, eam->rhor_spline);
|
||||||
interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline);
|
interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline);
|
||||||
|
|
||||||
// replicate data for multiple types;
|
// replicate data for multiple types;
|
||||||
for(int tt = 0; tt < ntypes * ntypes; tt++) {
|
for(int tt = 0; tt < ntypes * ntypes; tt++) {
|
||||||
for(int k = 0; k < nrho_tot; k++)
|
for(int k = 0; k < eam->nrho_tot; k++)
|
||||||
eam->frho_spline[tt*nrho_tot + k] = eam->frho_spline[k];
|
eam->frho_spline[tt*eam->nrho_tot + k] = eam->frho_spline[k];
|
||||||
for(int k = 0; k < nr_tot; k++)
|
for(int k = 0; k < eam->nr_tot; k++)
|
||||||
eam->rhor_spline[tt*nr_tot + k] = eam->rhor_spline[k];
|
eam->rhor_spline[tt*eam->nr_tot + k] = eam->rhor_spline[k];
|
||||||
for(int k = 0; k < nr_tot; k++)
|
for(int k = 0; k < eam->nr_tot; k++)
|
||||||
eam->z2r_spline[tt*nr_tot + k] = eam->z2r_spline[k];
|
eam->z2r_spline[tt*eam->nr_tot + k] = eam->z2r_spline[k];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) {
|
void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline) {
|
||||||
for(int m = 1; m <= n; m++) spline[m * 7 + 6] = f[m];
|
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[1 * 7 + 5] = spline[2 * 7 + 6] - spline[1 * 7 + 6];
|
||||||
@ -240,7 +247,7 @@ void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void grab(FILE* fptr, MMD_int n, MMD_float* list) {
|
void grab(FILE* fptr, int n, MD_FLOAT* list) {
|
||||||
char* ptr;
|
char* ptr;
|
||||||
char line[MAXLINE];
|
char line[MAXLINE];
|
||||||
int i = 0;
|
int i = 0;
|
||||||
|
@ -21,6 +21,7 @@
|
|||||||
* =======================================================================================
|
* =======================================================================================
|
||||||
*/
|
*/
|
||||||
#include <likwid-marker.h>
|
#include <likwid-marker.h>
|
||||||
|
#include <math.h>
|
||||||
|
|
||||||
#include <allocate.h>
|
#include <allocate.h>
|
||||||
#include <timing.h>
|
#include <timing.h>
|
||||||
@ -29,18 +30,21 @@
|
|||||||
#include <atom.h>
|
#include <atom.h>
|
||||||
#include <stats.h>
|
#include <stats.h>
|
||||||
#include <eam.h>
|
#include <eam.h>
|
||||||
|
#include <util.h>
|
||||||
|
|
||||||
double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) {
|
double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) {
|
||||||
if(atom->nmax > eam->nmax) {
|
if(eam->nmax < atom->Nmax) {
|
||||||
eam->nmax = atom->nmax;
|
eam->nmax = atom->Nmax;
|
||||||
if(eam->fp != NULL) { free(eam->fp); }
|
if(eam->fp != NULL) { free(eam->fp); }
|
||||||
eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->nmax * sizeof(MD_FLOAT));
|
eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT));
|
||||||
}
|
}
|
||||||
|
|
||||||
int Nlocal = atom->Nlocal;
|
int Nlocal = atom->Nlocal;
|
||||||
int* neighs;
|
int* neighs;
|
||||||
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT* fp = eam->fp;
|
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT* fp = eam->fp;
|
||||||
MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* fhro_spline = eam->fhro_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
|
MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
|
||||||
|
int rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; int rdrho = eam->rdrho;
|
||||||
|
int nrho = eam->nrho; int nrho_tot = eam->nrho_tot;
|
||||||
double S = getTimeStamp();
|
double S = getTimeStamp();
|
||||||
LIKWID_MARKER_START("force_eam_fp");
|
LIKWID_MARKER_START("force_eam_fp");
|
||||||
|
|
||||||
@ -67,7 +71,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
|
|||||||
|
|
||||||
if(rsq < cutforcesq) {
|
if(rsq < cutforcesq) {
|
||||||
MD_FLOAT p = sqrt(rsq) * rdr + 1.0;
|
MD_FLOAT p = sqrt(rsq) * rdr + 1.0;
|
||||||
int m = static_cast<int>(p);
|
int m = (int)(p);
|
||||||
m = m < nr - 1 ? m : nr - 1;
|
m = m < nr - 1 ? m : nr - 1;
|
||||||
p -= m;
|
p -= m;
|
||||||
p = p < 1.0 ? p : 1.0;
|
p = p < 1.0 ? p : 1.0;
|
||||||
@ -81,7 +85,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
|
|||||||
|
|
||||||
const int type_ii = type_i * type_i;
|
const int type_ii = type_i * type_i;
|
||||||
MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0;
|
MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0;
|
||||||
int m = static_cast<int>(p);
|
int m = (int)(p);
|
||||||
m = MAX(1, MIN(m, nrho - 1));
|
m = MAX(1, MIN(m, nrho - 1));
|
||||||
p -= m;
|
p -= m;
|
||||||
p = MIN(p, 1.0);
|
p = MIN(p, 1.0);
|
||||||
@ -116,9 +120,9 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
|
|||||||
const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
|
const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
|
||||||
|
|
||||||
if(rsq < cutforcesq) {
|
if(rsq < cutforcesq) {
|
||||||
MMD_float r = sqrt(rsq);
|
MD_FLOAT r = sqrt(rsq);
|
||||||
MMD_float p = r * rdr + 1.0;
|
MD_FLOAT p = r * rdr + 1.0;
|
||||||
int m = static_cast<int>(p);
|
int m = (int)(p);
|
||||||
m = m < nr - 1 ? m : nr - 1;
|
m = m < nr - 1 ? m : nr - 1;
|
||||||
p -= m;
|
p -= m;
|
||||||
p = p < 1.0 ? p : 1.0;
|
p = p < 1.0 ? p : 1.0;
|
||||||
@ -134,24 +138,24 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
|
|||||||
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
|
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
|
||||||
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
|
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
|
||||||
|
|
||||||
MMD_float rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p +
|
MD_FLOAT rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p +
|
||||||
rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p +
|
rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p +
|
||||||
rhor_spline[type_ij * nr_tot + m * 7 + 2];
|
rhor_spline[type_ij * nr_tot + m * 7 + 2];
|
||||||
|
|
||||||
MMD_float z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p +
|
MD_FLOAT z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p +
|
||||||
z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p +
|
z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p +
|
||||||
z2r_spline[type_ij * nr_tot + m * 7 + 2];
|
z2r_spline[type_ij * nr_tot + m * 7 + 2];
|
||||||
|
|
||||||
MMD_float z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p +
|
MD_FLOAT z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p +
|
||||||
z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p +
|
z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p +
|
||||||
z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p +
|
z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p +
|
||||||
z2r_spline[type_ij * nr_tot + m * 7 + 6];
|
z2r_spline[type_ij * nr_tot + m * 7 + 6];
|
||||||
|
|
||||||
MMD_float recip = 1.0 / r;
|
MD_FLOAT recip = 1.0 / r;
|
||||||
MMD_float phi = z2 * recip;
|
MD_FLOAT phi = z2 * recip;
|
||||||
MMD_float phip = z2p * recip - phi * recip;
|
MD_FLOAT phip = z2p * recip - phi * recip;
|
||||||
MMD_float psip = fp[i] * rhoip + fp[j] * rhoip + phip;
|
MD_FLOAT psip = fp[i] * rhoip + fp[j] * rhoip + phip;
|
||||||
MMD_float fpair = -psip * recip;
|
MD_FLOAT fpair = -psip * recip;
|
||||||
|
|
||||||
fix += delx * fpair;
|
fix += delx * fpair;
|
||||||
fiy += dely * fpair;
|
fiy += dely * fpair;
|
||||||
|
@ -30,14 +30,12 @@ typedef struct {
|
|||||||
MD_FLOAT *x, *y, *z;
|
MD_FLOAT *x, *y, *z;
|
||||||
MD_FLOAT *vx, *vy, *vz;
|
MD_FLOAT *vx, *vy, *vz;
|
||||||
MD_FLOAT *fx, *fy, *fz;
|
MD_FLOAT *fx, *fy, *fz;
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
int *type;
|
int *type;
|
||||||
int ntypes;
|
int ntypes;
|
||||||
MD_FLOAT *epsilon;
|
MD_FLOAT *epsilon;
|
||||||
MD_FLOAT *sigma6;
|
MD_FLOAT *sigma6;
|
||||||
MD_FLOAT *cutforcesq;
|
MD_FLOAT *cutforcesq;
|
||||||
MD_FLOAT *cutneighsq;
|
MD_FLOAT *cutneighsq;
|
||||||
#endif
|
|
||||||
} Atom;
|
} Atom;
|
||||||
|
|
||||||
extern void initAtom(Atom*);
|
extern void initAtom(Atom*);
|
||||||
|
@ -20,6 +20,8 @@
|
|||||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||||
* =======================================================================================
|
* =======================================================================================
|
||||||
*/
|
*/
|
||||||
|
#include <stdio.h>
|
||||||
|
|
||||||
#include <atom.h>
|
#include <atom.h>
|
||||||
#include <parameter.h>
|
#include <parameter.h>
|
||||||
|
|
||||||
@ -34,6 +36,7 @@ typedef struct {
|
|||||||
typedef struct {
|
typedef struct {
|
||||||
MD_FLOAT* fp;
|
MD_FLOAT* fp;
|
||||||
int nmax;
|
int nmax;
|
||||||
|
int ntypes;
|
||||||
int nrho, nr;
|
int nrho, nr;
|
||||||
int nrho_tot, nr_tot;
|
int nrho_tot, nr_tot;
|
||||||
MD_FLOAT dr, rdr, drho, rdrho;
|
MD_FLOAT dr, rdr, drho, rdrho;
|
||||||
@ -43,12 +46,12 @@ typedef struct {
|
|||||||
Funcfl* file;
|
Funcfl* file;
|
||||||
} Eam;
|
} Eam;
|
||||||
|
|
||||||
void init_eam(Eam* eam, int ntypes);
|
void initEam(Eam* eam, const char* input_file, int ntypes);
|
||||||
void coeff(Eam* eam, const char* arg);
|
void coeff(Eam* eam, const char* arg);
|
||||||
void init_style(Eam* eam);
|
void init_style(Eam* eam);
|
||||||
void read_file(Funcfl* file, const char* filename);
|
void read_file(Funcfl* file, const char* filename);
|
||||||
void file2array(Eam* eam);
|
void file2array(Eam* eam);
|
||||||
void array2spline(Eam* eam);
|
void array2spline(Eam* eam);
|
||||||
void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline);
|
void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline);
|
||||||
void grab(FILE* fptr, MMD_int n, MMD_float* list);
|
void grab(FILE* fptr, int n, MD_FLOAT* list);
|
||||||
#endif
|
#endif
|
||||||
|
@ -33,6 +33,7 @@
|
|||||||
#endif
|
#endif
|
||||||
|
|
||||||
typedef struct {
|
typedef struct {
|
||||||
|
int force_field;
|
||||||
char* input_file;
|
char* input_file;
|
||||||
MD_FLOAT epsilon;
|
MD_FLOAT epsilon;
|
||||||
MD_FLOAT sigma6;
|
MD_FLOAT sigma6;
|
||||||
|
@ -128,7 +128,6 @@ int main(int argc, const char *argv[]) {
|
|||||||
initAtom(atom);
|
initAtom(atom);
|
||||||
initStats(&stats);
|
initStats(&stats);
|
||||||
|
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->ntypes = param.ntypes;
|
atom->ntypes = param.ntypes;
|
||||||
atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
||||||
atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
|
||||||
@ -140,7 +139,6 @@ int main(int argc, const char *argv[]) {
|
|||||||
atom->cutneighsq[i] = param.cutneigh * param.cutneigh;
|
atom->cutneighsq[i] = param.cutneigh * param.cutneigh;
|
||||||
atom->cutforcesq[i] = param.cutforce * param.cutforce;
|
atom->cutforcesq[i] = param.cutforce * param.cutforce;
|
||||||
}
|
}
|
||||||
#endif
|
|
||||||
|
|
||||||
DEBUG("Creating atoms...\n");
|
DEBUG("Creating atoms...\n");
|
||||||
for(int i = 0; i < param.nx; ++i) {
|
for(int i = 0; i < param.nx; ++i) {
|
||||||
@ -176,9 +174,7 @@ int main(int argc, const char *argv[]) {
|
|||||||
for(int jj = 0; jj < fac_y; ++jj) {
|
for(int jj = 0; jj < fac_y; ++jj) {
|
||||||
for(int kk = 0; kk < fac_z; ++kk) {
|
for(int kk = 0; kk < fac_z; ++kk) {
|
||||||
if(added_atoms < atoms_per_unit_cell) {
|
if(added_atoms < atoms_per_unit_cell) {
|
||||||
#ifdef EXPLICIT_TYPES
|
|
||||||
atom->type[atom->Nlocal] = rand() % atom->ntypes;
|
atom->type[atom->Nlocal] = rand() % atom->ntypes;
|
||||||
#endif
|
|
||||||
ADD_ATOM(ii * offset_x, jj * offset_y, kk * offset_z, vx, vy, vz);
|
ADD_ATOM(ii * offset_x, jj * offset_y, kk * offset_z, vx, vy, vz);
|
||||||
added_atoms++;
|
added_atoms++;
|
||||||
}
|
}
|
||||||
|
Loading…
Reference in New Issue
Block a user