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

Addition of potential, code modifications to incorporate multicomponent spline MEAM in pair_meam_spline.

Browse Source 
Backwards compatible with previous version of pair_meam_spline.
This commit is contained in:
DallasTrinkle
2017-04-28 14:48:34 -05:00
parent 2f32fb7f8b
commit 7d9670bc6c
3 changed files with 566 additions and 211 deletions
Download Patch File Download Diff File Expand all files Collapse all files

130
potentials/TiO.meam.spline Normal file
View File

@ -0,0 +1,130 @@
# Ti-O cubic spline potential where O is in the dilute limit. DATE: 2016-06-05 CONTRIBUTOR: Pinchao Zhang, Dallas R. Trinkle
meam/spline 2 Ti O
spline3eq
13
-20 0
1.742692837 3.744277175966 99.4865081627958
2.05580176725 0.910839730906 10.8702523265355
2.3689106975 0.388045896634 -1.55322418749562
2.68201962775 -0.018840906533 2.43630041329215
2.995128558 -0.248098929639 2.67912713976835
3.30823748825 -0.264489550297 -0.125056384603077
3.6213464185 -0.227196189283 1.10662555360438
3.93445534875 -0.129293090176 -0.592053676745914
4.247564279 -0.059685366933 -0.470123414607672
4.56067320925 -0.031100025561 -0.0380739973059663
4.8737821395 -0.013847363202 -0.0711547960695406
5.18689106975 -0.003203412728 -0.081768292420175
5.5 0 -0.0571422964883619
spline3eq
5
0.155001355787331 0
1.9 0.533321679606674 0
2.8 0.456402081843862 -1.60311717015859
3.7 -0.324281383502201 1.19940299483249
4.6 -0.474029826906675 1.47909794595154
5.5 0 -2.49521499855605
spline3eq
13
0 0
1.742692837 0 0
2.05580176725 0 0
2.3689106975 0 0
2.68201962775 0 0
2.995128558 0 0
3.30823748825 0 0
3.6213464185 0 0
3.93445534875 0 0
4.247564279 0 0
4.56067320925 0 0
4.8737821395 0 0
5.18689106975 0 0
5.5 0 0
spline3eq
11
-1 0
2.055801767 1.7475279661 -525.869786904802
2.2912215903 -5.8677963945 252.796316927755
2.5266414136 -8.3376288737 71.7318388721015
2.7620612369 -5.8398712842 -1.93587742753693
2.9974810602 -3.1140648231 -39.2999192667503
3.2329008835 -1.7257245065 14.3424136002004
3.4683207068 -0.4428977017 -29.4925534559498
3.7037405301 -0.1466643003 -3.18010534572236
3.9391603534 -0.2095507945 3.33490838803603
4.1745801767 -0.1442384563 3.71918691359508
4.41 0 -9.66717019857564
spline3eq
5
-61.9827585211652 0
1.9 11.2293641315584 0
2.8 -27.9976343076148 122.648031332411
3.7 -8.32979773113248 -54.3340881766381
4.6 -1.00863195297399 3.23150064581724
5.5 0 -5.3514242228123
spline3eq
4
0.00776934946045395 0.105197706160344
-55.14233165 -0.29745568008 0.00152870603877451
-44.7409899033333 -0.15449458722 0.00038933722543571
-34.3396481566667 0.05098657168 0.00038124926922248
-23.93830641 0.57342694704 0.0156639264890892
spline3eq
5
-0.00676745157022662 -0.0159520381982146
-23.9928 0.297607384684645 0
-15.9241175 0.216691597077105 -0.0024248755353942
-7.855435 0.0637598673719069 0.00306245895013358
0.213247499999998 -0.00183450621970427 -0.00177588407633909
8.28193 -0.111277018874367 0
spline3eq
10
2.77327511656661 0
2.055801767 -0.1485215264 72.2010867146919
2.31737934844444 1.6845304918 -47.2744689053404
2.57895692988889 2.0113365977 -15.1859578405326
2.84053451133333 1.1444092747 3.33978204841873
3.10211209277778 0.2861606803 2.587867603808
3.36368967422222 -0.3459281126 6.14070694084556
3.62526725566667 -0.6257480601 3.7397696717154
3.88684483711111 -0.6119510826 4.64749084871402
4.14842241855556 -0.3112059651 2.83275746415936
4.41 0 -15.0612086827734
spline3eq
5
12.3315547862781 0
1.9 2.62105440156724 0
2.8 10.2850803058354 -25.439802988016
3.7 3.23933763743897 -7.20203673434025
4.6 -5.79049355858613 39.5509978688682
5.5 0 -41.221771373642
spline3eq
8
8.33642274810572 -60.4024574736564
-1 0.07651409193 -110.652321293778
-0.724509054371429 0.14155824541 44.8853405500508
-0.449018108742857 0.75788697341 -25.3065115342002
-0.173527163114286 0.63011570378 -2.48510144915082
0.101963782514286 0.09049597305 2.68769386908235
0.377454728142857 -0.35741586657 -1.01558570129633
0.652945673771428 -0.65293217647 13.4224786001212
0.9284366194 -6.00912190653 -452.752542694929
spline3eq
5
0.137191606537625 -1.55094230968985
-1 0.0513843442016519 0
-0.5 0.0179024412245673 -2.44986494990154
0 -0.260650876879273 3.91774583656401
0.5 -0.190163791764901 -4.84414871911743
1 -0.763795416646599 0
spline3eq
8
0 0
-1 0 0
-0.724509054371429 0 0
-0.449018108742857 0 0
-0.173527163114286 0 0
0.101963782514286 0 0
0.377454728142857 0 0
0.652945673771428 0 0
0.9284366194 0 0

517
src/USER-MISC/pair_meam_spline.cpp
View File

@ -13,6 +13,8 @@
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Contributing author: Alexander Stukowski (LLNL), alex@stukowski.com Contributing author: Alexander Stukowski (LLNL), alex@stukowski.com
Will Tipton (Cornell), wwt26@cornell.edu
Dallas R. Trinkle (UIUC), dtrinkle@illinois.edu / Pinchao Zhang (UIUC)
see LLNL copyright notice at bottom of file see LLNL copyright notice at bottom of file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -23,12 +25,14 @@
* 25-Mar-11 - AS: Fixed calculation of per-atom virial stress. * 25-Mar-11 - AS: Fixed calculation of per-atom virial stress.
* 11-Apr-11 - AS: Adapted code to new memory management of LAMMPS. * 11-Apr-11 - AS: Adapted code to new memory management of LAMMPS.
* 24-Sep-11 - AS: Adapted code to new interface of Error::one() function. * 24-Sep-11 - AS: Adapted code to new interface of Error::one() function.
* 20-Jun-13 - WT: Added support for multiple species types
* 25-Apr-17 - DRT/PZ: Modified format of multiple species type to conform with pairing
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include <math.h> #include "math.h"
#include <stdio.h> #include "stdio.h"
#include <stdlib.h> #include "stdlib.h"
#include <string.h> #include "string.h"
#include "pair_meam_spline.h" #include "pair_meam_spline.h"
#include "atom.h" #include "atom.h"
#include "force.h" #include "force.h"
@ -39,6 +43,9 @@
#include "neigh_request.h" #include "neigh_request.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include <iostream>
using namespace std;
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -49,7 +56,6 @@ PairMEAMSpline::PairMEAMSpline(LAMMPS *lmp) : Pair(lmp)
single_enable = 0; single_enable = 0;
restartinfo = 0; restartinfo = 0;
one_coeff = 1; one_coeff = 1;
manybody_flag = 1;
nelements = 0; nelements = 0;
elements = NULL; elements = NULL;
@ -77,6 +83,15 @@ PairMEAMSpline::~PairMEAMSpline()
if(allocated) { if(allocated) {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(cutsq); memory->destroy(cutsq);
delete[] phis;
delete[] Us;
delete[] rhos;
delete[] fs;
delete[] gs;
delete[] zero_atom_energies;
delete [] map; delete [] map;
} }
} }
@ -85,11 +100,12 @@ PairMEAMSpline::~PairMEAMSpline()
void PairMEAMSpline::compute(int eflag, int vflag) void PairMEAMSpline::compute(int eflag, int vflag)
{ {
if (eflag || vflag) ev_setup(eflag, vflag); if (eflag || vflag) {
else evflag = vflag_fdotr = ev_setup(eflag, vflag);
eflag_global = vflag_global = eflag_atom = vflag_atom = 0; } else {
evflag = vflag_fdotr = eflag_global = 0;
double cutforcesq = cutoff*cutoff; vflag_global = eflag_atom = vflag_atom = 0;
}
// Grow per-atom array if necessary // Grow per-atom array if necessary
@ -99,22 +115,13 @@ void PairMEAMSpline::compute(int eflag, int vflag)
memory->create(Uprime_values,nmax,"pair:Uprime"); memory->create(Uprime_values,nmax,"pair:Uprime");
} }
double** const x = atom->x;
double** forces = atom->f;
int nlocal = atom->nlocal;
bool newton_pair = force->newton_pair;
int inum_full = listfull->inum;
int* ilist_full = listfull->ilist;
int* numneigh_full = listfull->numneigh;
int** firstneigh_full = listfull->firstneigh;
// Determine the maximum number of neighbors a single atom has // Determine the maximum number of neighbors a single atom has
int newMaxNeighbors = 0; int newMaxNeighbors = 0;
for(int ii = 0; ii < inum_full; ii++) { for(int ii = 0; ii < listfull->inum; ii++) {
int jnum = numneigh_full[ilist_full[ii]]; int jnum = listfull->numneigh[listfull->ilist[ii]];
if(jnum > newMaxNeighbors) newMaxNeighbors = jnum; if(jnum > newMaxNeighbors)
newMaxNeighbors = jnum;
} }
// Allocate array for temporary bond info // Allocate array for temporary bond info
@ -126,35 +133,88 @@ void PairMEAMSpline::compute(int eflag, int vflag)
} }
// Sum three-body contributions to charge density and // Sum three-body contributions to charge density and
// compute embedding energies // the embedding energy
for(int ii = 0; ii < inum_full; ii++) { for(int ii = 0; ii < listfull->inum; ii++) {
int i = ilist_full[ii]; int i = listfull->ilist[ii];
double xtmp = x[i][0];
double ytmp = x[i][1];
double ztmp = x[i][2];
int* jlist = firstneigh_full[i];
int jnum = numneigh_full[i];
double rho_value = 0;
int numBonds = 0; int numBonds = 0;
// compute charge density and numBonds
double rho_value = compute_three_body_contrib_to_charge_density(i, numBonds);
//cout<<"i "<<i<<" rho "<<rho_value<<" numBonds "<<numBonds<<endl;
// Compute embedding energy and its derivative
double Uprime_i = compute_embedding_energy_and_deriv(eflag, i, rho_value);
//cout<<"i "<<i<<" U prime "<<Uprime_i<<endl;
// Compute three-body contributions to force
compute_three_body_contrib_to_forces(i, numBonds, Uprime_i);
}
// Communicate U'(rho) values
comm->forward_comm_pair(this);
// Compute two-body pair interactions
compute_two_body_pair_interactions();
if(vflag_fdotr)
virial_fdotr_compute();
}
double PairMEAMSpline::pair_density(int i)
{
double rho_value = 0;
MEAM2Body* nextTwoBodyInfo = twoBodyInfo; MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
for(int jj = 0; jj < jnum; jj++) { for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
int j = jlist[jj]; int j = listfull->firstneigh[i][jj];
j &= NEIGHMASK; j &= NEIGHMASK;
double jdelx = x[j][0] - xtmp; double jdelx = atom->x[j][0] - atom->x[i][0];
double jdely = x[j][1] - ytmp; double jdely = atom->x[j][1] - atom->x[i][1];
double jdelz = x[j][2] - ztmp; double jdelz = atom->x[j][2] - atom->x[i][2];
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
double rij = sqrt(rij_sq);
if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq);
rho_value += rhos[i_to_potl(j)].eval(rij);
}
}
return rho_value;
}
double PairMEAMSpline::three_body_density(int i)
{
double rho_value = 0;
int numBonds=0;
MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
int j = listfull->firstneigh[i][jj];
j &= NEIGHMASK;
double jdelx = atom->x[j][0] - atom->x[i][0];
double jdely = atom->x[j][1] - atom->x[i][1];
double jdelz = atom->x[j][2] - atom->x[i][2];
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz; double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
if(rij_sq < cutforcesq) { if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq); double rij = sqrt(rij_sq);
double partial_sum = 0; double partial_sum = 0;
nextTwoBodyInfo->tag = j; nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij; nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime); nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij; nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij; nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij; nextTwoBodyInfo->del[2] = jdelz / rij;
@ -164,31 +224,82 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] + double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] + nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]); nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * g.eval(cos_theta); partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
} }
rho_value += nextTwoBodyInfo->f * partial_sum; rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rho.eval(rij); numBonds++;
nextTwoBodyInfo++;
}
}
return rho_value;
}
double PairMEAMSpline::compute_three_body_contrib_to_charge_density(int i, int& numBonds) {
double rho_value = 0;
MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
int j = listfull->firstneigh[i][jj];
j &= NEIGHMASK;
double jdelx = atom->x[j][0] - atom->x[i][0];
double jdely = atom->x[j][1] - atom->x[i][1];
double jdelz = atom->x[j][2] - atom->x[i][2];
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq);
double partial_sum = 0;
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = twoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rhos[i_to_potl(j)].eval(rij);
numBonds++; numBonds++;
nextTwoBodyInfo++; nextTwoBodyInfo++;
} }
} }
// Compute embedding energy and its derivative return rho_value;
}
double PairMEAMSpline::compute_embedding_energy_and_deriv(int eflag, int i, double rho_value) {
double Uprime_i; double Uprime_i;
double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy; double embeddingEnergy = Us[i_to_potl(i)].eval(rho_value, Uprime_i)
- zero_atom_energies[i_to_potl(i)];
//cout<<"Density "<<rho_value<<endl;
//cout<<"embedding energy "<<embeddingEnergy<<endl;
Uprime_values[i] = Uprime_i; Uprime_values[i] = Uprime_i;
if(eflag) { if(eflag) {
if(eflag_global) eng_vdwl += embeddingEnergy; if(eflag_global)
if(eflag_atom) eatom[i] += embeddingEnergy; eng_vdwl += embeddingEnergy;
if(eflag_atom)
eatom[i] += embeddingEnergy;
} }
return Uprime_i;
}
void PairMEAMSpline::compute_three_body_contrib_to_forces(int i, int numBonds, double Uprime_i) {
double forces_i[3] = {0, 0, 0}; double forces_i[3] = {0, 0, 0};
// Compute three-body contributions to force
for(int jj = 0; jj < numBonds; jj++) { for(int jj = 0; jj < numBonds; jj++) {
const MEAM2Body bondj = twoBodyInfo[jj]; const MEAM2Body bondj = twoBodyInfo[jj];
double rij = bondj.r; double rij = bondj.r;
@ -207,7 +318,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
bondj.del[1]*bondk->del[1] + bondj.del[1]*bondk->del[1] +
bondj.del[2]*bondk->del[2]); bondj.del[2]*bondk->del[2]);
double g_prime; double g_prime;
double g_value = g.eval(cos_theta, g_prime); double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
double f_rik_prime = bondk->fprime; double f_rik_prime = bondk->fprime;
double f_rik = bondk->f; double f_rik = bondk->f;
@ -237,9 +348,9 @@ void PairMEAMSpline::compute(int eflag, int vflag)
forces_i[2] -= fk[2]; forces_i[2] -= fk[2];
int k = bondk->tag; int k = bondk->tag;
forces[k][0] += fk[0]; atom->f[k][0] += fk[0];
forces[k][1] += fk[1]; atom->f[k][1] += fk[1];
forces[k][2] += fk[2]; atom->f[k][2] += fk[2];
if(evflag) { if(evflag) {
double delta_ij[3]; double delta_ij[3];
@ -254,76 +365,91 @@ void PairMEAMSpline::compute(int eflag, int vflag)
} }
} }
forces[i][0] -= forces_j[0]; atom->f[i][0] -= forces_j[0];
forces[i][1] -= forces_j[1]; atom->f[i][1] -= forces_j[1];
forces[i][2] -= forces_j[2]; atom->f[i][2] -= forces_j[2];
forces[j][0] += forces_j[0]; atom->f[j][0] += forces_j[0];
forces[j][1] += forces_j[1]; atom->f[j][1] += forces_j[1];
forces[j][2] += forces_j[2]; atom->f[j][2] += forces_j[2];
} }
forces[i][0] += forces_i[0]; atom->f[i][0] += forces_i[0];
forces[i][1] += forces_i[1]; atom->f[i][1] += forces_i[1];
forces[i][2] += forces_i[2]; atom->f[i][2] += forces_i[2];
} }
// Communicate U'(rho) values void PairMEAMSpline::compute_two_body_pair_interactions() {
for(int ii = 0; ii < listhalf->inum; ii++) {
int i = listhalf->ilist[ii];
comm->forward_comm_pair(this); for(int jj = 0; jj < listhalf->numneigh[i]; jj++) {
int j = listhalf->firstneigh[i][jj];
int inum_half = listhalf->inum;
int* ilist_half = listhalf->ilist;
int* numneigh_half = listhalf->numneigh;
int** firstneigh_half = listhalf->firstneigh;
// Compute two-body pair interactions
for(int ii = 0; ii < inum_half; ii++) {
int i = ilist_half[ii];
double xtmp = x[i][0];
double ytmp = x[i][1];
double ztmp = x[i][2];
int* jlist = firstneigh_half[i];
int jnum = numneigh_half[i];
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK; j &= NEIGHMASK;
double jdel[3]; double jdel[3];
jdel[0] = x[j][0] - xtmp; jdel[0] = atom->x[j][0] - atom->x[i][0];
jdel[1] = x[j][1] - ytmp; jdel[1] = atom->x[j][1] - atom->x[i][1];
jdel[2] = x[j][2] - ztmp; jdel[2] = atom->x[j][2] - atom->x[i][2];
double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2]; double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2];
if(rij_sq < cutforcesq) { if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq); double rij = sqrt(rij_sq);
double rho_prime; // double rho_prime;
rho.eval(rij, rho_prime); // rhos[i_to_potl(j)].eval(rij, rho_prime);
double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]); // cout<<"rho_prime "<<rho_prime<<endl;
// double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]);
double rho_prime_i,rho_prime_j;
rhos[i_to_potl(i)].eval(rij,rho_prime_i);
rhos[i_to_potl(j)].eval(rij,rho_prime_j);
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
//cout<<"fpair "<<fpair<<endl;
double pair_pot_deriv; double pair_pot_deriv;
double pair_pot = phi.eval(rij, pair_pot_deriv); double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
//cout<<"pair potential "<<pair_pot<<endl;
//cout<<"pair pot_deriv "<<pair_pot_deriv<<endl;
fpair += pair_pot_deriv; fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient // Divide by r_ij to get forces from gradient
fpair /= rij; fpair /= rij;
forces[i][0] += jdel[0]*fpair; atom->f[i][0] += jdel[0]*fpair;
forces[i][1] += jdel[1]*fpair; //cout<<i<<"pair fx "<<atom->f[i][0]<<endl;
forces[i][2] += jdel[2]*fpair; atom->f[i][1] += jdel[1]*fpair;
forces[j][0] -= jdel[0]*fpair; //cout<<i<<"pair fy "<<atom->f[i][1]<<endl;
forces[j][1] -= jdel[1]*fpair; atom->f[i][2] += jdel[2]*fpair;
forces[j][2] -= jdel[2]*fpair; //cout<<i<<"pair fz "<<atom->f[i][2]<<endl;
if (evflag) ev_tally(i, j, nlocal, newton_pair, atom->f[j][0] -= jdel[0]*fpair;
//cout<<j<<"pair fx "<<atom->f[j][0]<<endl;
atom->f[j][1] -= jdel[1]*fpair;
//cout<<j<<"pair fy "<<atom->f[j][1]<<endl;
atom->f[j][2] -= jdel[2]*fpair;
////cout<<j<<"pair fz "<<atom->f[j][2]<<endl;
if (evflag) ev_tally(i, j, atom->nlocal, force->newton_pair,
pair_pot, 0.0, -fpair, jdel[0], jdel[1], jdel[2]); pair_pot, 0.0, -fpair, jdel[0], jdel[1], jdel[2]);
} }
} }
} }
}
if(vflag_fdotr) virial_fdotr_compute(); /* ----------------------------------------------------------------------
helper functions to map atom types to potential array indices
------------------------------------------------------------------------- */
int PairMEAMSpline::ij_to_potl(int i, int j) {
int n = atom->ntypes;
int itype = atom->type[i];
int jtype = atom->type[j];
// printf("%d %d %d\n",n,itype,jtype);
return jtype - 1 + (itype-1)*n - (itype-1)*itype/2;
}
int PairMEAMSpline::i_to_potl(int i) {
int itype = atom->type[i];
return itype - 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -331,11 +457,23 @@ void PairMEAMSpline::compute(int eflag, int vflag)
void PairMEAMSpline::allocate() void PairMEAMSpline::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ntypes; int n = nelements;
memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cutsq,n+1,n+1,"pair:cutsq");
int nmultichoose2 = n*(n+1)/2;
//Change the functional form
//f_ij->f_i
//g_i(cos\theta_ijk)->g_jk(cos\theta_ijk)
phis = new SplineFunction[nmultichoose2];
Us = new SplineFunction[n];
rhos = new SplineFunction[n];
fs = new SplineFunction[n];
gs = new SplineFunction[nmultichoose2];
zero_atom_energies = new double[n];
map = new int[n+1]; map = new int[n+1];
} }
@ -356,7 +494,8 @@ void PairMEAMSpline::coeff(int narg, char **arg)
{ {
int i,j,n; int i,j,n;
if (!allocated) allocate(); if (!allocated)
allocate();
if (narg != 3 + atom->ntypes) if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients"); error->all(FLERR,"Incorrect args for pair coefficients");
@ -366,45 +505,34 @@ void PairMEAMSpline::coeff(int narg, char **arg)
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients"); error->all(FLERR,"Incorrect args for pair coefficients");
// read potential file: also sets the number of elements.
read_file(arg[2]);
// read args that map atom types to elements in potential file // read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL // map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements // nelements = # of unique elements
// elements = list of element names // elements = list of element names
if (elements) { if ((nelements == 1) && (strlen(elements[0]) == 0)) {
for (i = 0; i < nelements; i++) delete [] elements[i]; // old style: we only have one species, so we're either "NULL" or we match.
delete [] elements; for (i = 3; i < narg; i++)
} if (strcmp(arg[i],"NULL") == 0)
elements = new char*[atom->ntypes]; map[i-2] = -1;
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL; else
map[i-2] = 0;
nelements = 0; } else {
for (i = 3; i < narg; i++) { for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) { if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1; map[i-2] = -1;
continue; continue;
} }
for (j = 0; j < nelements; j++) for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break; if (strcmp(arg[i],elements[j]) == 0)
map[i-2] = j; break;
if (j == nelements) { if (j < nelements) map[i-2] = j;
n = strlen(arg[i]) + 1; else error->all(FLERR,"No matching element in EAM potential file");
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
} }
} }
// for now, only allow single element
if (nelements > 1)
error->all(FLERR,
"Pair meam/spline only supports single element potentials");
// read potential file
read_file(arg[2]);
// clear setflag since coeff() called once with I,J = * * // clear setflag since coeff() called once with I,J = * *
n = atom->ntypes; n = atom->ntypes;
@ -425,57 +553,126 @@ void PairMEAMSpline::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
} }
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
#define MAXLINE 1024 #define MAXLINE 1024
void PairMEAMSpline::read_file(const char* filename) void PairMEAMSpline::read_file(const char* filename)
{ {
int nmultichoose2; // = (n+1)*n/2;
if(comm->me == 0) { if(comm->me == 0) {
FILE *fp = force->open_potential(filename); FILE *fp = fopen(filename, "r");
if(fp == NULL) { if(fp == NULL) {
char str[1024]; char str[1024];
sprintf(str,"Cannot open spline MEAM potential file %s", filename); sprintf(str,"Cannot open spline MEAM potential file %s", filename);
error->one(FLERR,str); error->one(FLERR,str);
} }
// Skip first line of file. // Skip first line of file. It's a comment.
char line[MAXLINE]; char line[MAXLINE];
fgets(line, MAXLINE, fp); fgets(line, MAXLINE, fp);
// Second line holds potential type (currently just "meam/spline") in new potential format.
bool isNewFormat;
long loc = ftell(fp);
fgets(line, MAXLINE, fp);
if (strncmp(line, "meam/spline", 11) == 0) {
isNewFormat = true;
// parse the rest of the line!
char *linep = line+12, *word;
const char *sep = " ,;:-\t\n"; // overkill, but safe
word = strsep(&linep, sep);
if (! *word)
error->one(FLERR, "Need to include number of atomic species on meam/spline line in potential file");
int n = atoi(word);
if (n<1)
error->one(FLERR, "Invalid number of atomic species on meam/spline line in potential file");
nelements = n;
elements = new char*[n];
for (int i=0; i<n; ++i) {
word = strsep(&linep, sep);
if (! *word)
error->one(FLERR, "Not enough atomic species in meam/spline\n");
elements[i] = new char[strlen(word)+1];
strcpy(elements[i], word);
}
} else {
isNewFormat = false;
nelements = 1; // old format only handles one species anyway; this is for backwards compatibility
elements = new char*[1];
elements[0] = new char[1];
strcpy(elements[0], "");
fseek(fp, loc, SEEK_SET);
}
nmultichoose2 = ((nelements+1)*nelements)/2;
// allocate!!
allocate();
// Parse spline functions. // Parse spline functions.
phi.parse(fp, error);
rho.parse(fp, error); for (int i = 0; i < nmultichoose2; i++)
U.parse(fp, error); phis[i].parse(fp, error, isNewFormat);
f.parse(fp, error); for (int i = 0; i < nelements; i++)
g.parse(fp, error); rhos[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nelements; i++)
Us[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nelements; i++)
fs[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nmultichoose2; i++)
gs[i].parse(fp, error, isNewFormat);
fclose(fp); fclose(fp);
} }
// Transfer spline functions from master processor to all other processors. // Transfer spline functions from master processor to all other processors.
phi.communicate(world, comm->me); MPI_Bcast(&nelements, 1, MPI_INT, 0, world);
rho.communicate(world, comm->me); MPI_Bcast(&nmultichoose2, 1, MPI_INT, 0, world);
f.communicate(world, comm->me); // allocate!!
U.communicate(world, comm->me); if (!allocated) {
g.communicate(world, comm->me); allocate();
elements = new char*[nelements];
}
for (int i = 0; i < nelements; ++i) {
int n;
if (comm->me == 0)
n = strlen(elements[i]);
MPI_Bcast(&n, 1, MPI_INT, 0, world);
if (comm->me != 0)
elements[i] = new char[n];
MPI_Bcast(elements[i], n, MPI_CHAR, 0, world);
}
for (int i = 0; i < nmultichoose2; i++)
phis[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
rhos[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
fs[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
Us[i].communicate(world, comm->me);
for (int i = 0; i < nmultichoose2; i++)
gs[i].communicate(world, comm->me);
// Calculate 'zero-point energy' of single atom in vacuum. // Calculate 'zero-point energy' of single atom in vacuum.
zero_atom_energy = U.eval(0.0); for (int i = 0; i < nelements; i++)
zero_atom_energies[i] = Us[i].eval(0.0);
// Determine maximum cutoff radius of all relevant spline functions. // Determine maximum cutoff radius of all relevant spline functions.
cutoff = 0.0; cutoff = 0.0;
if(phi.cutoff() > cutoff) cutoff = phi.cutoff(); for (int i = 0; i < nmultichoose2; i++)
if(rho.cutoff() > cutoff) cutoff = rho.cutoff(); if(phis[i].cutoff() > cutoff)
if(f.cutoff() > cutoff) cutoff = f.cutoff(); cutoff = phis[i].cutoff();
for (int i = 0; i < nelements; i++)
if(rhos[i].cutoff() > cutoff)
cutoff = rhos[i].cutoff();
for (int i = 0; i < nelements; i++)
if(fs[i].cutoff() > cutoff)
cutoff = fs[i].cutoff();
// Set LAMMPS pair interaction flags. // Set LAMMPS pair interaction flags.
for(int i = 1; i <= atom->ntypes; i++) { for(int i = 1; i <= atom->ntypes; i++) {
for(int j = 1; j <= atom->ntypes; j++) { for(int j = 1; j <= atom->ntypes; j++) {
setflag[i][j] = 1; // setflag[i][j] = 1;
cutsq[i][j] = cutoff; cutsq[i][j] = cutoff; // should this be squared?
} }
} }
@ -495,12 +692,15 @@ void PairMEAMSpline::init_style()
error->all(FLERR,"Pair style meam/spline requires newton pair on"); error->all(FLERR,"Pair style meam/spline requires newton pair on");
// Need both full and half neighbor list. // Need both full and half neighbor list.
int irequest_full = neighbor->request(this,instance_me); int irequest_full = neighbor->request(this);
neighbor->requests[irequest_full]->id = 1; neighbor->requests[irequest_full]->id = 1;
neighbor->requests[irequest_full]->half = 0; neighbor->requests[irequest_full]->half = 0;
neighbor->requests[irequest_full]->full = 1; neighbor->requests[irequest_full]->full = 1;
int irequest_half = neighbor->request(this,instance_me); int irequest_half = neighbor->request(this);
neighbor->requests[irequest_half]->id = 2; neighbor->requests[irequest_half]->id = 2;
// neighbor->requests[irequest_half]->half = 0;
// neighbor->requests[irequest_half]->half_from_full = 1;
// neighbor->requests[irequest_half]->otherlist = irequest_full;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -523,14 +723,13 @@ double PairMEAMSpline::init_one(int i, int j)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int PairMEAMSpline::pack_forward_comm(int n, int *list, double *buf, int PairMEAMSpline::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{ {
int* list_iter = list; int* list_iter = list;
int* list_iter_end = list + n; int* list_iter_end = list + n;
while(list_iter != list_iter_end) while(list_iter != list_iter_end)
*buf++ = Uprime_values[*list_iter++]; *buf++ = Uprime_values[*list_iter++];
return n; return 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -563,10 +762,14 @@ double PairMEAMSpline::memory_usage()
/// Parses the spline knots from a text file. /// Parses the spline knots from a text file.
void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error) void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error, bool isNewFormat)
{ {
char line[MAXLINE]; char line[MAXLINE];
// If new format, read the spline format. Should always be "spline3eq" for now.
if (isNewFormat)
fgets(line, MAXLINE, fp);
// Parse number of spline knots. // Parse number of spline knots.
fgets(line, MAXLINE, fp); fgets(line, MAXLINE, fp);
int n = atoi(line); int n = atoi(line);
@ -578,8 +781,10 @@ void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error)
double d0 = atof(strtok(line, " \t\n\r\f")); double d0 = atof(strtok(line, " \t\n\r\f"));
double dN = atof(strtok(NULL, " \t\n\r\f")); double dN = atof(strtok(NULL, " \t\n\r\f"));
init(n, d0, dN); init(n, d0, dN);
//printf("%s\n",line);
// Skip line. // Skip line in old format
if (!isNewFormat)
fgets(line, MAXLINE, fp); fgets(line, MAXLINE, fp);
// Parse knot coordinates. // Parse knot coordinates.
@ -734,3 +939,5 @@ void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename, const ch
* Lawrence Livermore National Security, LLC, and shall not be used for * Lawrence Livermore National Security, LLC, and shall not be used for
* advertising or product endorsement purposes. * advertising or product endorsement purposes.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

50
src/USER-MISC/pair_meam_spline.h
View File

@ -1,4 +1,4 @@
/* -*- c++ -*- ---------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
@ -43,10 +43,24 @@ public:
virtual void compute(int, int); virtual void compute(int, int);
void settings(int, char **); void settings(int, char **);
void coeff(int, char **); void coeff(int, char **);
void get_coeff(double *, double *);
double pair_density(int );
double three_body_density(int );
void init_style(); void init_style();
void init_list(int, class NeighList *); void init_list(int, class NeighList *);
double init_one(int, int); double init_one(int, int);
// helper functions for compute()
double compute_three_body_contrib_to_charge_density(int i, int& numBonds); // returns rho_value and returns numBonds by reference
double compute_embedding_energy_and_deriv(int eflag, int i, double rho_value); // returns the derivative of the embedding energy Uprime_i
void compute_three_body_contrib_to_forces(int i, int numBonds, double Uprime_i);
void compute_two_body_pair_interactions();
int ij_to_potl(int i, int j);
int i_to_potl(int i);
int pack_forward_comm(int, int *, double *, int, int *); int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *); void unpack_forward_comm(int, int, double *);
int pack_reverse_comm(int, int, double *); int pack_reverse_comm(int, int, double *);
@ -74,15 +88,15 @@ protected:
} }
/// Initialization of spline function. /// Initialization of spline function.
void init(int _n, double _deriv0, double _derivN) { void init(int _N, double _deriv0, double _derivN) {
N = _n; N = _N;
deriv0 = _deriv0; deriv0 = _deriv0;
derivN = _derivN; derivN = _derivN;
delete[] X; // if (X) delete[] X;
delete[] Xs; // if (Xs) delete[] Xs;
delete[] Y; // if (Y) delete[] Y;
delete[] Y2; // if (Y2) delete[] Y2;
delete[] Ydelta; // if (Ydelta) delete[] Ydelta;
X = new double[N]; X = new double[N];
Xs = new double[N]; Xs = new double[N];
Y = new double[N]; Y = new double[N];
@ -97,7 +111,7 @@ protected:
int numKnots() const { return N; } int numKnots() const { return N; }
/// Parses the spline knots from a text file. /// Parses the spline knots from a text file.
void parse(FILE* fp, Error* error); void parse(FILE* fp, Error* error, bool isNewFormat);
/// Calculates the second derivatives of the cubic spline. /// Calculates the second derivatives of the cubic spline.
void prepareSpline(Error* error); void prepareSpline(Error* error);
@ -209,18 +223,18 @@ protected:
/// Helper data structure for potential routine. /// Helper data structure for potential routine.
struct MEAM2Body { struct MEAM2Body {
int tag; int tag; // holds the index of the second atom (j)
double r; double r;
double f, fprime; double f, fprime;
double del[3]; double del[3];
}; };
SplineFunction phi; // Phi(r_ij) SplineFunction* phis; // Phi_i(r_ij)
SplineFunction rho; // Rho(r_ij) SplineFunction* rhos; // Rho_ij(r_ij)
SplineFunction f; // f(r_ij) SplineFunction* fs; // f_i(r_ij)
SplineFunction U; // U(rho) SplineFunction* Us; // U_i(rho)
SplineFunction g; // g(cos_theta) SplineFunction* gs; // g_ij(cos_theta)
double zero_atom_energy; // Shift embedding energy by this value to make it zero for a single atom in vacuum. double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
double cutoff; // The cutoff radius double cutoff; // The cutoff radius
@ -231,6 +245,8 @@ protected:
void read_file(const char* filename); void read_file(const char* filename);
void allocate(); void allocate();
}; };
} }
@ -279,3 +295,5 @@ protected:
* *
* See file 'pair_spline_meam.cpp' for history of changes. * See file 'pair_spline_meam.cpp' for history of changes.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */