Convert to C++, allow multiple instances (03c93b31) · Commits · 郑智淋 / lammps

src/USER-MEAMC/README

+1 −1

Original line number	Diff line number	Diff line
		@@ -2,7 +2,7 @@ This package implements the MEAM/C potential as a LAMMPS pair style.

		+==============================================================================+

		This package is a translation of the MEAM package to native C. See
		This package is a translation of the MEAM package to native C++. See
		that package as well as the Fortran code distributed in lib/meam for
		the original sources and information.

src/USER-MEAMC/meam.h

+186 −124

Original line number	Diff line number	Diff line
		#ifndef LMP_MEAM_H
		#define LMP_MEAM_H

		#include <stdlib.h>

		#define maxelt 5

		extern "C" {
		double fm_exp(double);
		}

		typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;

		typedef struct {
		typedef struct
		{
		int dim1, dim2;
		double* ptr;
		} allocatable_double_2d;

		typedef struct {
		class MEAM {
		private:
		// cutforce = force cutoff
		// cutforcesq = force cutoff squared

		@@ -58,7 +66,8 @@ typedef struct {

		double Ec_meam[maxelt + 1][maxelt + 1], re_meam[maxelt + 1][maxelt + 1];
		double Omega_meam[maxelt + 1], Z_meam[maxelt + 1];
		double A_meam[maxelt+1],alpha_meam[maxelt+1][maxelt+1],rho0_meam[maxelt+1];
		double A_meam[maxelt + 1], alpha_meam[maxelt + 1][maxelt + 1],
		rho0_meam[maxelt + 1];
		double delta_meam[maxelt + 1][maxelt + 1];
		double beta0_meam[maxelt + 1], beta1_meam[maxelt + 1];
		double beta2_meam[maxelt + 1], beta3_meam[maxelt + 1];
		@@ -74,9 +83,11 @@ typedef struct {

		allocatable_double_2d phir; // [:,:]

		allocatable_double_2d phirar,phirar1,phirar2,phirar3,phirar4,phirar5,phirar6; // [:,:]
		allocatable_double_2d phirar, phirar1, phirar2, phirar3, phirar4, phirar5,
		phirar6; // [:,:]

		double attrac_meam[maxelt+1][maxelt+1],repuls_meam[maxelt+1][maxelt+1];
		double attrac_meam[maxelt + 1][maxelt + 1],
		repuls_meam[maxelt + 1][maxelt + 1];

		double Cmin_meam[maxelt + 1][maxelt + 1][maxelt + 1];
		double Cmax_meam[maxelt + 1][maxelt + 1][maxelt + 1];
		@@ -89,29 +100,77 @@ typedef struct {

		int nr, nrar;
		double dr, rdrar;
		} meam_data_t;

		extern meam_data_t meam_data;
		public:
		void meam_checkindex(int, int, int, int, int);
		void meam_setup_param_(int, double, int, int, int*);
		void meam_dens_final_(int, int, int, int, int, double, double, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, int);
		void G_gam(double, int, double, double, int);
		void dG_gam(double, int, double, double, double);
		void meam_dens_init_(int, int, int, int, int, double, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, int);
		void getscreen(int, int, double, double, double, double, int, int, int, int, int, int, int);
		void screen(int, int, int, double, double, double, int, int, int, int, int*);
		void calc_rho1(int, int, int, int, int, double, int, int, double, double, double, double, double, double, double, double, double, double);
		void dsij(int, int, int, int, int, int, double, double, double, int, int, int, double, double, double*);
		void fcut(double, double*);
		void dfcut(double, double, double);
		void dCfunc(double, double, double, double*);
		void dCfunc2(double, double, double, double, double);

		void meam_force_(int, int, int, int, int, int, double, double, int, int, int, double, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, int*);
		void meam_cleanup_();
		void meam_setup_done_(double*);
		void alloyparams();
		void compute_pair_meam();
		double phi_meam(double, int, int);
		void compute_reference_density();
		void get_shpfcn(double*, lattice_t);
		void get_tavref(double, double, double, double, double, double, double, double, double, double, double, double, double, int, int, lattice_t);
		void get_Zij(int*, lattice_t);
		void get_Zij2(int, double, double*, lattice_t, double, double);
		void get_sijk(double, int, int, int, double*);
		void get_densref(double, int, int, double, double, double, double, double, double, double, double);
		double zbl(double, int, int);
		double erose(double, double, double, double, double, double, int);
		void interpolate_meam(int);
		double compute_phi(double, int, int);
		void meam_setup_global_(int, int, double, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, int*);

		};

		// Functions we need for compat
		#ifndef max
		#define max(a, b) \
		({ __typeof__ (a) _a = (a); \
		({ \
		__typeof__(a) _a = (a); \
		__typeof__(b) _b = (b); \
		_a > _b ? _a : _b; })
		_a > _b ? _a : _b; \
		})
		#endif

		#ifndef min
		#define min(a, b) \
		({ __typeof__ (a) _a = (a); \
		({ \
		__typeof__(a) _a = (a); \
		__typeof__(b) _b = (b); \
		_a < _b ? _a : _b; })
		_a < _b ? _a : _b; \
		})
		#endif

		#define iszero(f) (fabs(f) < 1e-20)

		#define setall2d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) arr[__i][__j] = v;}
		#define setall3d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) for(int __k=1;__k<=maxelt;__k++) arr[__i][__j][__k] = v;}
		#define setall2d(arr, v) \
		{ \
		for (int __i = 1; __i <= maxelt; __i++) \
		for (int __j = 1; __j <= maxelt; __j++) \
		arr[__i][__j] = v; \
		}
		#define setall3d(arr, v) \
		{ \
		for (int __i = 1; __i <= maxelt; __i++) \
		for (int __j = 1; __j <= maxelt; __j++) \
		for (int __k = 1; __k <= maxelt; __k++) \
		arr[__i][__j][__k] = v; \
		}

		/*
		Fortran Array Semantics in C.
		@@ -125,33 +184,36 @@ Fortran Array Semantics in C.
		// we receive a pointer to the first element, and F dimensions is ptr(a,b,c)
		// we know c data structure is ptr[c][b][a]
		#define arrdim2v(ptr, a, b) \
		const int DIM1__##ptr = a; const int DIM2__##ptr = b;\
		(void)(DIM1__##ptr); (void)(DIM2__##ptr);
		const int DIM1__##ptr = a; \
		const int DIM2__##ptr = b; \
		(void)(DIM1__##ptr); \
		(void)(DIM2__##ptr);
		#define arrdim3v(ptr, a, b, c) \
		const int DIM1__##ptr = a; const int DIM2__##ptr = b; const int DIM3__##ptr = c;\
		const int DIM1__##ptr = a; \
		const int DIM2__##ptr = b; \
		const int DIM3__##ptr = c; \
		(void)(DIM1__##ptr); (void)(DIM2__##ptr; (void)(DIM3__##ptr);


		// access data with same index as used in fortran (1-based)
		#define arr1v(ptr,i) \
		ptr[i-1]
		#define arr2v(ptr,i,j) \
		ptr[(DIM1__##ptr)*(j-1) + (i-1)]
		#define arr1v(ptr, i) ptr[i - 1]
		#define arr2v(ptr, i, j) ptr[(DIM1__##ptr) * (j - 1) + (i - 1)]
		#define arr3v(ptr, i, j, k) \
		ptr[(i-1) + (j-1)(DIM1__##ptr) + (k-1)(DIM1__##ptr)*(DIM2__##ptr)]
		ptr[(i - 1) + (j - 1) * (DIM1__##ptr) + \
		(k - 1) * (DIM1__##ptr) * (DIM2__##ptr)]

		// allocatable arrays via special type
		#define allocate_2d(arr, cols, rows) \
		arr.dim1 = cols; arr.dim2=rows; arr.ptr=(double)calloc((size_t)(cols)(size_t)(rows),sizeof(double));
		#define allocated(a) \
		(a.ptr!=NULL)
		#define deallocate(a) \
		({free(a.ptr); a.ptr=NULL;})
		arr.dim1 = cols; \
		arr.dim2 = rows; \
		arr.ptr = (double)calloc((size_t)(cols) (size_t)(rows), sizeof(double));
		#define allocated(a) (a.ptr != NULL)
		#define meam_deallocate(a) \
		({ \
		free(a.ptr); \
		a.ptr = NULL; \
		})
		// access data with same index as used in fortran (1-based)
		#define arr2(arr,i,j) \
		arr.ptr[(arr.dim1)*(j-1) + (i-1)]



		#define arr2(arr, i, j) arr.ptr[(arr.dim1) * (j - 1) + (i - 1)]


		#endif
		No newline at end of file

src/USER-MEAMC/meam_cleanup.cpp

+9 −11

Original line number	Diff line number	Diff line
		extern "C" {
		#include "meam.h"

		void
		meam_cleanup_(void)
		MEAM::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);
		}
		meam_deallocate(this->phirar6);
		meam_deallocate(this->phirar5);
		meam_deallocate(this->phirar4);
		meam_deallocate(this->phirar3);
		meam_deallocate(this->phirar2);
		meam_deallocate(this->phirar1);
		meam_deallocate(this->phirar);
		meam_deallocate(this->phir);
		}
		No newline at end of file

src/USER-MEAMC/meam_data.cpp

deleted100644 → 0

+0 −5

Original line number	Diff line number	Diff line
		extern "C" {
		#include "meam.h"

		meam_data_t meam_data = {};
		}
		No newline at end of file

src/USER-MEAMC/meam_dens_final.cpp

+35 −45

Original line number	Diff line number	Diff line
		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);

		// in meam_setup_done
		void get_shpfcn(double* s /* s(3) */, lattice_t latt);

		// Extern "C" declaration has the form:
		//
		// void meam_dens_final_(int , int , int , int , int , double , double *,
		@@ -28,7 +19,7 @@ void get_shpfcn(double* s /* s(3) */, lattice_t latt);
		//

		void
		meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		MEAM::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,
		@@ -67,23 +58,23 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		for (m = 1; m <= 6; m++) {
		arr1v(rho2, i) =
		arr1v(rho2, i) +
		meam_data.v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i);
		this->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);
		this->v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i);
		}

		if (arr1v(rho0, i) > 0.0) {
		if (meam_data.ialloy == 1) {
		if (this->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 if (this->ialloy == 2) {
		arr2v(t_ave, 1, i) = this->t1_meam[elti];
		arr2v(t_ave, 2, i) = this->t2_meam[elti];
		arr2v(t_ave, 3, i) = this->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);
		@@ -99,56 +90,56 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i));
		}

		Z = meam_data.Z_meam[elti];
		Z = this->Z_meam[elti];

		G_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti],
		meam_data.gsmooth_factor, &G, errorflag);
		G_gam(arr1v(Gamma, i), this->ibar_meam[elti],
		this->gsmooth_factor, &G, errorflag);
		if (*errorflag != 0)
		return;
		get_shpfcn(shp, meam_data.lattce_meam[elti][elti]);
		if (meam_data.ibar_meam[elti] <= 0) {
		get_shpfcn(shp, this->lattce_meam[elti][elti]);
		if (this->ibar_meam[elti] <= 0) {
		Gbar = 1.0;
		dGbar = 0.0;
		} else {
		if (meam_data.mix_ref_t == 1) {
		if (this->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]) /
		gam = (this->t1_meam[elti] * shp[1] +
		this->t2_meam[elti] * shp[2] +
		this->t3_meam[elti] * shp[3]) /
		(Z * Z);
		}
		G_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor, &Gbar,
		G_gam(gam, this->ibar_meam[elti], this->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) {
		if (this->mix_ref_t == 1) {
		if (this->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,
		dG_gam(gam, this->ibar_meam[elti], this->gsmooth_factor,
		&Gbar, &dGbar);
		}
		rho_bkgd = meam_data.rho0_meam[elti] * Z * Gbar;
		rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
		} else {
		if (meam_data.bkgd_dyn == 1) {
		rho_bkgd = meam_data.rho0_meam[elti] * Z;
		if (this->bkgd_dyn == 1) {
		rho_bkgd = this->rho0_meam[elti] * Z;
		} else {
		rho_bkgd = meam_data.rho_ref_meam[elti];
		rho_bkgd = this->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);
		dG_gam(arr1v(Gamma, i), this->ibar_meam[elti],
		this->gsmooth_factor, &G, &dG);

		arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom;

		@@ -161,30 +152,30 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		// 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) {
		if (this->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];
		B = this->A_meam[elti] * this->Ec_meam[elti][elti];

		if (!iszero(rhob)) {
		if (meam_data.emb_lin_neg == 1 && rhob <= 0) {
		if (this->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) {
		if (this->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) {
		if (this->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);
		@@ -192,7 +183,7 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		}
		}
		} else {
		if (meam_data.emb_lin_neg == 1) {
		if (this->emb_lin_neg == 1) {
		arr1v(fp, i) = -B;
		} else {
		arr1v(fp, i) = B;
		@@ -205,7 +196,7 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
		// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

		void
		G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag)
		MEAM::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)
		@@ -246,7 +237,7 @@ G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag)
		// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

		void
		dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG)
		MEAM::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)
		@@ -289,4 +280,3 @@ dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG)
		}

		// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
		}
		No newline at end of file

Admin message