Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65680352
pair_lcbop.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jun 5, 12:12
Size
40 KB
Mime Type
text/x-c
Expires
Fri, Jun 7, 12:12 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18100634
Attached To
rLAMMPS lammps
pair_lcbop.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Dominik Wójt (Wroclaw University of Technology)
based on pair_airebo by Ase Henry (MIT)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "mpi.h"
#include "pair_lcbop.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "my_page.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define TOL 1.0e-9
#define PGDELTA 1
/* ---------------------------------------------------------------------- */
PairLCBOP::PairLCBOP(LAMMPS *lmp) : Pair(lmp) {
single_enable = 0;
one_coeff = 1;
manybody_flag = 1;
ghostneigh = 1;
maxlocal = 0;
SR_numneigh = NULL;
SR_firstneigh = NULL;
ipage = NULL;
pgsize = oneatom = 0;
N = NULL;
M = NULL;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairLCBOP::~PairLCBOP()
{
memory->destroy(SR_numneigh);
memory->sfree(SR_firstneigh);
delete [] ipage;
memory->destroy(N);
memory->destroy(M);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cutghost);
delete [] map;
}
}
/* ---------------------------------------------------------------------- */
void PairLCBOP::compute(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = vflag_atom = 0;
SR_neigh();
FSR(eflag,vflag);
FLR(eflag,vflag);
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLCBOP::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutghost,n+1,n+1,"pair:cutghost");
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLCBOP::settings(int narg, char **arg) {
if( narg != 0 ) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLCBOP::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to C and NULL
// map[i] = which element (0 for C) the Ith atom type is, -1 if NULL
for (int i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
} else if (strcmp(arg[i],"C") == 0) {
map[i-2] = 0;
} else error->all(FLERR,"Incorrect args for pair coefficients");
}
// read potential file and initialize fitting splines
read_file(arg[2]);
spline_init();
// clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLCBOP::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style LCBOP requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style LCBOP requires newton pair on");
// need a full neighbor list, including neighbors of ghosts
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->ghost = 1;
// local SR neighbor list
// create pages if first time or if neighbor pgsize/oneatom has changed
int create = 0;
if (ipage == NULL) create = 1;
if (pgsize != neighbor->pgsize) create = 1;
if (oneatom != neighbor->oneatom) create = 1;
if (create) {
delete [] ipage;
pgsize = neighbor->pgsize;
oneatom = neighbor->oneatom;
int nmypage = comm->nthreads;
ipage = new MyPage<int>[nmypage];
for (int i = 0; i < nmypage; i++)
ipage[i].init(oneatom,pgsize,PGDELTA);
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLCBOP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
// cut3rebo = 3 SR distances
cut3rebo = 3.0 * r_2;
// cutmax = furthest distance from an owned atom
// at which another atom will feel force, i.e. the ghost cutoff
// for SR term in potential:
// interaction = M-K-I-J-L-N with I = owned and J = ghost
// I to N is max distance = 3 SR distances
// for V_LR term in potential:
// r_2_LR
// cutghost = SR cutoff used in SR_neigh() for neighbors of ghosts
double cutmax = MAX( cut3rebo,r_2_LR );
cutghost[i][j] = r_2;
cutLRsq = r_2_LR*r_2_LR;
cutghost[j][i] = cutghost[i][j];
r_2_sq = r_2*r_2;
return cutmax;
}
/* ----------------------------------------------------------------------
create SR neighbor list from main neighbor list
SR neighbor list stores neighbors of ghost atoms
------------------------------------------------------------------------- */
void PairLCBOP::SR_neigh()
{
int i,j,ii,jj,n,allnum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS;
int *ilist,*jlist,*numneigh,**firstneigh;
int *neighptr;
double **x = atom->x;
if (atom->nmax > maxlocal) { // ensure ther is enough space
maxlocal = atom->nmax; // for atoms and ghosts allocated
memory->destroy(SR_numneigh);
memory->sfree(SR_firstneigh);
memory->destroy(N);
memory->destroy(M);
memory->create(SR_numneigh,maxlocal,"LCBOP:numneigh");
SR_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),
"LCBOP:firstneigh");
memory->create(N,maxlocal,"LCBOP:N");
memory->create(M,maxlocal,"LCBOP:M");
}
allnum = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// store all SR neighs of owned and ghost atoms
// scan full neighbor list of I
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
n = 0;
neighptr = ipage->vget();
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
N[i] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < r_2_sq) {
neighptr[n++] = j;
N[i] += f_c(sqrt(rsq),r_1,r_2,&dS);
}
}
SR_firstneigh[i] = neighptr;
SR_numneigh[i] = n;
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
// calculate M_i
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
M[i] = 0.0;
jlist = SR_firstneigh[i];
jnum = SR_numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < r_2_sq) {
double f_c_ij = f_c(sqrt(rsq),r_1,r_2,&dS);
double Nji = N[j]-f_c_ij;
// F(xij) = 1-f_c_LR(Nji, 2,3,&dummy)
M[i] += f_c_ij * ( 1-f_c_LR(Nji, 2,3,&dS) );
}
}
}
}
/* ----------------------------------------------------------------------
Short range forces and energy
------------------------------------------------------------------------- */
void PairLCBOP::FSR(int eflag, int vflag)
{
int i,j,jj,ii,inum;
tagint itag,jtag;
double delx,dely,delz,fpair,xtmp,ytmp,ztmp;
double r_sq,rijmag,f_c_ij,df_c_ij;
double VR,dVRdi,VA,Bij,dVAdi,dVA;
double del[3];
int *ilist,*SR_neighs;
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
// two-body interactions from SR neighbor list, skip half of them
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
SR_neighs = SR_firstneigh[i];
for (jj = 0; jj < SR_numneigh[i]; jj++) {
j = SR_neighs[jj];
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
r_sq = delx*delx + dely*dely + delz*delz;
rijmag = sqrt(r_sq);
f_c_ij = f_c( rijmag,r_1,r_2,&df_c_ij );
if( f_c_ij <= TOL ) continue;
VR = A*exp(-alpha*rijmag);
dVRdi = -alpha*VR;
dVRdi = dVRdi*f_c_ij + df_c_ij*VR; // VR -> VR * f_c_ij
VR *= f_c_ij;
VA = dVA = 0.0;
{
double term = B_1 * exp(-beta_1*rijmag);
VA += term;
dVA += -beta_1 * term;
term = B_2 * exp(-beta_2*rijmag);
VA += term;
dVA += -beta_2 * term;
}
dVA = dVA*f_c_ij + df_c_ij*VA; // VA -> VA * f_c_ij
VA *= f_c_ij;
del[0] = delx;
del[1] = dely;
del[2] = delz;
Bij = bondorder(i,j,del,rijmag,VA,f,vflag_atom);
dVAdi = Bij*dVA;
// F = (dVRdi+dVAdi)*(-grad rijmag)
// grad_i rijmag = \vec{rij} /rijmag
// grad_j rijmag = -\vec{rij} /rijmag
fpair = -(dVRdi-dVAdi) / rijmag;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
double evdwl=0.0;
if (eflag) evdwl = VR - Bij*VA;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
/* ----------------------------------------------------------------------
compute long range forces and energy
------------------------------------------------------------------------- */
void PairLCBOP::FLR(int eflag, int vflag)
{
int i,j,jj,ii;
tagint itag,jtag;
double delx,dely,delz,fpair,xtmp,ytmp,ztmp;
double r_sq,rijmag,f_c_ij,df_c_ij;
double V,dVdi;
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int inum = list->inum;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
// two-body interactions from full neighbor list, skip half of them
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
int *neighs = firstneigh[i];
for (jj = 0; jj < numneigh[i]; jj++) {
j = neighs[jj];
j &= NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
r_sq = delx*delx + dely*dely + delz*delz;
rijmag = sqrt(r_sq);
f_c_ij = 1-f_c( rijmag,r_1,r_2,&df_c_ij );
df_c_ij = -df_c_ij;
// derivative may be inherited from previous call, see f_c_LR definition
f_c_ij *= f_c_LR( rijmag, r_1_LR, r_2_LR, &df_c_ij );
if( f_c_ij <= TOL ) continue;
V = dVdi = 0;
if( rijmag<r_0 ) {
double exp_part = exp( -lambda_1*(rijmag-r_0) );
V = eps_1*( exp_part*exp_part - 2*exp_part) + v_1;
dVdi = 2*eps_1*lambda_1*exp_part*( 1-exp_part );
} else {
double exp_part = exp( -lambda_2*(rijmag-r_0) );
V = eps_2*( exp_part*exp_part - 2*exp_part) + v_2;
dVdi = 2*eps_2*lambda_2*exp_part*( 1-exp_part );
}
dVdi = dVdi*f_c_ij + df_c_ij*V; // V -> V * f_c_ij
V *= f_c_ij;
// F = (dVdi)*(-grad rijmag)
// grad_i rijmag = \vec{rij} /rijmag
// grad_j rijmag = -\vec{rij} /rijmag
fpair = -dVdi / rijmag;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
double evdwl=0.0;
if (eflag) evdwl = V;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
/* ----------------------------------------------------------------------
forces for Nij and Mij
------------------------------------------------------------------------- */
void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom ) {
int atomi = i;
int atomj = j;
int *SR_neighs = SR_firstneigh[i];
double **x = atom->x;
for( int k=0; k<SR_numneigh[i]; k++ ) {
int atomk = SR_neighs[k];
if (atomk != atomj) {
double rik[3];
rik[0] = x[atomi][0]-x[atomk][0];
rik[1] = x[atomi][1]-x[atomk][1];
rik[2] = x[atomi][2]-x[atomk][2];
double riksq = (rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2]);
if( riksq > r_1*r_1 ) { // && riksq < r_2*r_2, if second condition not fulfilled neighbor would not be in the list
double rikmag = sqrt(riksq);
double df_c_ik;
f_c( rikmag, r_1, r_2, &df_c_ik );
// F = factor*df_c_ik*(-grad rikmag)
// grad_i rikmag = \vec{rik} /rikmag
// grad_k rikmag = -\vec{rik} /rikmag
double fpair = -factor*df_c_ik / rikmag;
f[atomi][0] += rik[0]*fpair;
f[atomi][1] += rik[1]*fpair;
f[atomi][2] += rik[2]*fpair;
f[atomk][0] -= rik[0]*fpair;
f[atomk][1] -= rik[1]*fpair;
f[atomk][2] -= rik[2]*fpair;
if (vflag_atom) v_tally2(atomi,atomk,fpair,rik);
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom ) {
int atomi = i;
int atomj = j;
int *SR_neighs = SR_firstneigh[i];
double **x = atom->x;
for( int k=0; k<SR_numneigh[i]; k++ ) {
int atomk = SR_neighs[k];
if (atomk != atomj) {
double rik[3];
rik[0] = x[atomi][0]-x[atomk][0];
rik[1] = x[atomi][1]-x[atomk][1];
rik[2] = x[atomi][2]-x[atomk][2];
double rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2]));
double df_c_ik;
double f_c_ik = f_c( rikmag, r_1, r_2, &df_c_ik );
double Nki = N[k]-(f_c_ik);
// double Mij = M[i] - f_c_ij*( 1-f_c(Nji, 2,3,&dummy) );
double dF=0;
double Fx = 1-f_c_LR(Nki, 2,3,&dF);
dF = -dF;
if( df_c_ik > TOL ) {
double factor2 = factor*df_c_ik*Fx;
// F = factor2*(-grad rikmag)
// grad_i rikmag = \vec{rik} /rikmag
// grad_k rikmag = -\vec{rik} /rikmag
double fpair = -factor2 / rikmag;
f[atomi][0] += rik[0]*fpair;
f[atomi][1] += rik[1]*fpair;
f[atomi][2] += rik[2]*fpair;
f[atomk][0] -= rik[0]*fpair;
f[atomk][1] -= rik[1]*fpair;
f[atomk][2] -= rik[2]*fpair;
if (vflag_atom) v_tally2(atomi,atomk,fpair,rik);
}
if( dF > TOL ) {
double factor2 = factor*f_c_ik*dF;
FNij( atomk, atomi, factor2, f, vflag_atom );
}
}
}
}
/* ----------------------------------------------------------------------
Bij function
------------------------------------------------------------------------- */
double PairLCBOP::bondorder(int i, int j, double rij[3],
double rijmag, double VA,
double **f, int vflag_atom)
{
double bij, bji;
/* bij & bji */{
double rji[3];
rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
bij = b(i,j,rij,rijmag,VA,f,vflag_atom);
bji = b(j,i,rji,rijmag,VA,f,vflag_atom);
}
double Fij_conj;
/* F_conj */{
double dummy;
double df_c_ij;
double f_c_ij = f_c( rijmag, r_1, r_2, &df_c_ij );
double Nij = MIN( 3, N[i]-(f_c_ij) );
double Nji = MIN( 3, N[j]-(f_c_ij) );
// F(xij) = 1-f_c(Nji, 2,3,&dummy)
double Mij = M[i] - f_c_ij*( 1-f_c(Nji, 2,3,&dummy) );
double Mji = M[j] - f_c_ij*( 1-f_c(Nij, 2,3,&dummy) );
Mij = MIN( Mij, 3 );
Mji = MIN( Mji, 3 );
double Nij_el, dNij_el_dNij, dNij_el_dMij;
double Nji_el, dNji_el_dNji, dNji_el_dMji;
{
double num_Nij_el = 4 - Mij;
double num_Nji_el = 4 - Mji;
double den_Nij_el = Nij + 1 - Mij;
double den_Nji_el = Nji + 1 - Mji;
Nij_el = num_Nij_el / den_Nij_el;
Nji_el = num_Nji_el / den_Nji_el;
dNij_el_dNij = -Nij_el/den_Nij_el;
dNji_el_dNji = -Nji_el/den_Nji_el;
dNij_el_dMij = ( -1 + Nij_el ) /den_Nij_el;
dNji_el_dMji = ( -1 + Nji_el ) /den_Nji_el;
}
double Nconj;
double dNconj_dNij;
double dNconj_dNji;
double dNconj_dNel;
{
double num_Nconj = ( Nij+1 )*( Nji+1 )*( Nij_el+Nji_el ) - 4*( Nij+Nji+2);
double den_Nconj = Nij*( 3-Nij )*( Nji+1 ) + Nji*( 3-Nji )*( Nij+1 ) + eps;
Nconj = num_Nconj / den_Nconj;
if( Nconj <= 0 ) {
Nconj = 0;
dNconj_dNij = 0;
dNconj_dNji = 0;
dNconj_dNel = 0;
} else if( Nconj >= 1 ) {
Nconj = 1;
dNconj_dNij = 0;
dNconj_dNji = 0;
dNconj_dNel = 0;
} else {
dNconj_dNij = (
( (Nji+1)*(Nij_el + Nji_el)-4)
- Nconj*( (Nji+1)*(3-2*Nij) + Nji*(3-Nji) )
) /den_Nconj;
dNconj_dNji = (
( (Nij+1)*(Nji_el + Nij_el)-4)
- Nconj*( (Nij+1)*(3-2*Nji) + Nij*(3-Nij) )
) /den_Nconj;
dNconj_dNel = (Nij+1)*(Nji+1) / den_Nconj;
}
}
double dF_dNij, dF_dNji, dF_dNconj;
Fij_conj = F_conj( Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj );
/*forces for Nij*/
if( 3-Nij > TOL ) {
double factor = -VA*0.5*( dF_dNij + dF_dNconj*( dNconj_dNij + dNconj_dNel*dNij_el_dNij ) );
FNij( i, j, factor, f, vflag_atom );
}
/*forces for Nji*/
if( 3-Nji > TOL ) {
double factor = -VA*0.5*( dF_dNji + dF_dNconj*( dNconj_dNji + dNconj_dNel*dNji_el_dNji ) );
FNij( j, i, factor, f, vflag_atom );
}
/*forces for Mij*/
if( 3-Mij > TOL ) {
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNij_el_dMij );
FMij( i, j, factor, f, vflag_atom );
}
if( 3-Mji > TOL ) {
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNji_el_dMji );
FMij( j, i, factor, f, vflag_atom );
}
}
double Bij = 0.5*( bij + bji + Fij_conj );
return Bij;
}
/* ----------------------------------------------------------------------
bij function
------------------------------------------------------------------------- */
double PairLCBOP::b(int i, int j, double rij[3],
double rijmag, double VA,
double **f, int vflag_atom) {
int *SR_neighs = SR_firstneigh[i];
double **x = atom->x;
int atomi = i;
int atomj = j;
//calculate bij magnitude
double bij = 1.0;
for (int k = 0; k < SR_numneigh[i]; k++) {
int atomk = SR_neighs[k];
if (atomk != atomj) {
double rik[3];
rik[0] = x[atomi][0]-x[atomk][0];
rik[1] = x[atomi][1]-x[atomk][1];
rik[2] = x[atomi][2]-x[atomk][2];
double rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2]));
double delta_ijk = rijmag-rikmag;
double dummy;
double f_c_ik = f_c( rikmag, r_1, r_2, &dummy );
double cos_ijk = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2]))
/ (rijmag*rikmag);
cos_ijk = MIN(cos_ijk,1.0);
cos_ijk = MAX(cos_ijk,-1.0);
double G = gSpline(cos_ijk, &dummy);
double H = hSpline(delta_ijk, &dummy);
bij += (f_c_ik*G*H);
}
}
bij = pow( bij, -delta );
// bij forces
for (int k = 0; k < SR_numneigh[i]; k++) {
int atomk = SR_neighs[k];
if (atomk != atomj) {
double rik[3];
rik[0] = x[atomi][0]-x[atomk][0];
rik[1] = x[atomi][1]-x[atomk][1];
rik[2] = x[atomi][2]-x[atomk][2];
double rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2]));
double delta_ijk = rijmag-rikmag;
double df_c_ik;
double f_c_ik = f_c( rikmag, r_1, r_2, &df_c_ik );
double cos_ijk = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2]))
/ (rijmag*rikmag);
cos_ijk = MIN(cos_ijk,1.0);
cos_ijk = MAX(cos_ijk,-1.0);
double dcos_ijk_dri[3],dcos_ijk_drj[3],dcos_ijk_drk[3];
dcos_ijk_drj[0] = -rik[0] / (rijmag*rikmag)
+ cos_ijk * rij[0] / (rijmag*rijmag);
dcos_ijk_drj[1] = -rik[1] / (rijmag*rikmag)
+ cos_ijk * rij[1] / (rijmag*rijmag);
dcos_ijk_drj[2] = -rik[2] / (rijmag*rikmag)
+ cos_ijk * rij[2] / (rijmag*rijmag);
dcos_ijk_drk[0] = -rij[0] / (rijmag*rikmag)
+ cos_ijk * rik[0] / (rikmag*rikmag);
dcos_ijk_drk[1] = -rij[1] / (rijmag*rikmag)
+ cos_ijk * rik[1] / (rikmag*rikmag);
dcos_ijk_drk[2] = -rij[2] / (rijmag*rikmag)
+ cos_ijk * rik[2] / (rikmag*rikmag);
dcos_ijk_dri[0] = -dcos_ijk_drk[0] - dcos_ijk_drj[0];
dcos_ijk_dri[1] = -dcos_ijk_drk[1] - dcos_ijk_drj[1];
dcos_ijk_dri[2] = -dcos_ijk_drk[2] - dcos_ijk_drj[2];
double dG, dH;
double G = gSpline( cos_ijk, &dG );
double H = hSpline( delta_ijk, &dH );
double tmp = -VA*0.5*(-0.5*bij*bij*bij);
double fi[3], fj[3], fk[3];
double tmp2 = -tmp*df_c_ik*G*H/rikmag;
// F = tmp*df_c_ik*G*H*(-grad rikmag)
// grad_i rikmag = \vec{rik} /rikmag
// grad_k rikmag = -\vec{rik} /rikmag
fi[0] = tmp2*rik[0];
fi[1] = tmp2*rik[1];
fi[2] = tmp2*rik[2];
fk[0] = -tmp2*rik[0];
fk[1] = -tmp2*rik[1];
fk[2] = -tmp2*rik[2];
tmp2 = -tmp*f_c_ik*dG*H;
// F = tmp*f_c_ik*dG*H*(-grad cos_ijk)
// grad_i cos_ijk = dcos_ijk_dri
// grad_j cos_ijk = dcos_ijk_drj
// grad_k cos_ijk = dcos_ijk_drk
fi[0] += tmp2*dcos_ijk_dri[0];
fi[1] += tmp2*dcos_ijk_dri[1];
fi[2] += tmp2*dcos_ijk_dri[2];
fj[0] = tmp2*dcos_ijk_drj[0];
fj[1] = tmp2*dcos_ijk_drj[1];
fj[2] = tmp2*dcos_ijk_drj[2];
fk[0] += tmp2*dcos_ijk_drk[0];
fk[1] += tmp2*dcos_ijk_drk[1];
fk[2] += tmp2*dcos_ijk_drk[2];
tmp2 = -tmp*f_c_ik*G*dH;
// F = tmp*f_c_ik*G*dH*(-grad delta_ijk)
// grad_i delta_ijk = \vec{rij} /rijmag - \vec{rik} /rijmag
// grad_j delta_ijk = -\vec{rij} /rijmag
// grad_k delta_ijk = \vec{rik} /rikmag
fi[0] += tmp2*( rij[0]/rijmag - rik[0]/rikmag );
fi[1] += tmp2*( rij[1]/rijmag - rik[1]/rikmag );
fi[2] += tmp2*( rij[2]/rijmag - rik[2]/rikmag );
fj[0] += tmp2*( -rij[0]/rijmag );
fj[1] += tmp2*( -rij[1]/rijmag );
fj[2] += tmp2*( -rij[2]/rijmag );
fk[0] += tmp2*( rik[0]/rikmag );
fk[1] += tmp2*( rik[1]/rikmag );
fk[2] += tmp2*( rik[2]/rikmag );
f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2];
f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2];
if (vflag_atom) {
double rji[3], rki[3];
rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];
v_tally3(atomi,atomj,atomk,fj,fk,rji,rki);
}
}
}
return bij;
}
/* ----------------------------------------------------------------------
spline interpolation for G
------------------------------------------------------------------------- */
void PairLCBOP::g_decompose_x( double x, size_t *field_idx, double *offset ) {
size_t i=0;
while( i<(6-1) && !( x<gX[i+1] ) )
i++;
*field_idx = i;
*offset = ( x - gX[i] );
}
/* ---------------------------------------------------------------------- */
double PairLCBOP::gSpline( double x, double *dgdc ) {
size_t i;
double x_n;
g_decompose_x( x, &i, &x_n );
double sum = 0;
*dgdc = 0;
double pow_x_n = 1.0;
for( size_t j=0; j<5; j++ ) {
sum += gC[j][i]*pow_x_n;
*dgdc += gC[j+1][i]*(j+1)*pow_x_n;
pow_x_n *= x_n;
}
sum += gC[5][i]*pow_x_n;
return sum;
}
/* ---------------------------------------------------------------------- */
double PairLCBOP::hSpline( double x, double *dhdx ) {
if( x < -d ) {
double z = kappa*( x+d );
double y = pow(z, 10.0);
double w = pow( 1+y, -0.1 );
*dhdx = kappa*L*w/(1+y);
return L*( 1 + z*w );
}
if( x > d ) {
*dhdx = R_1;
return R_0 + R_1*( x-d );
}
double result = 1 + C_1*x;
*dhdx = C_1*result;
double pow_x = x*x;
result += 0.5*C_1*C_1*pow_x;
pow_x *= x;// == x^3
*dhdx += 4*C_4*pow_x;
pow_x *= x;// == x^4
result += C_4*pow_x;
pow_x *= x;// == x^5
*dhdx += 6*C_6*pow_x;
pow_x *= x;// == x^5
result += C_6*pow_x;
return result;
}
/* ---------------------------------------------------------------------- */
double PairLCBOP::F_conj( double N_ij, double N_ji, double N_conj_ij, double *dFN_ij, double *dFN_ji, double *dFN_ij_conj ) {
size_t N_ij_int = MIN( static_cast<size_t>( floor( N_ij ) ), 2 ); // 2 is the highest number of field
size_t N_ji_int = MIN( static_cast<size_t>( floor( N_ji ) ), 2 ); // cast to suppress warning
double x = N_ij - N_ij_int;
double y = N_ji - N_ji_int;
const TF_conj_field &f0 = F_conj_field[N_ij_int][N_ji_int][0];
const TF_conj_field &f1 = F_conj_field[N_ij_int][N_ji_int][1];
double F_0 = 0;
double F_1 = 0;
double dF_0_dx = 0, dF_0_dy = 0;
double dF_1_dx = 0, dF_1_dy = 0;
double l, r;
if( N_conj_ij < 1 ) {
l = (1-y)* (1-x); r = ( f0.f_00 + x* x* f0.f_x_10 + y* y* f0.f_y_01 ); F_0 += l*r; dF_0_dx += -(1-y)*r +l*2*x* f0.f_x_10; dF_0_dy += -(1-x)*r +l*2*y* f0.f_y_01;
l = (1-y)* x; r = ( f0.f_10 + (1-x)*(1-x)*f0.f_x_00 + y* y* f0.f_y_11 ); F_0 += l*r; dF_0_dx += (1-y)*r -l*2*(1-x)*f0.f_x_00; dF_0_dy += -x* r +l*2*y* f0.f_y_11;
l = y* (1-x); r = ( f0.f_01 + x* x* f0.f_x_11 + (1-y)*(1-y)*f0.f_y_00 ); F_0 += l*r; dF_0_dx += -y* r +l*2*x* f0.f_x_11; dF_0_dy += (1-x)*r -l*2*(1-y)*f0.f_y_00;
l = y* x; r = ( f0.f_11 + (1-x)*(1-x)*f0.f_x_01 + (1-y)*(1-y)*f0.f_y_10 ); F_0 += l*r; dF_0_dx += y* r -l*2*(1-x)*f0.f_x_01; dF_0_dy += x* r -l*2*(1-y)*f0.f_y_10;
}
if( N_conj_ij > 0 ) {
l = (1-y)* (1-x); r = ( f0.f_00 + x* x* f1.f_x_10 + y* y* f1.f_y_01 ); F_1 += l*r; dF_1_dx += -(1-y)*r +l*2*x* f1.f_x_10; dF_1_dy += -(1-x)*r +l*2*y* f1.f_y_01;
l = (1-y)* x; r = ( f1.f_10 + (1-x)*(1-x)*f1.f_x_00 + y* y* f1.f_y_11 ); F_1 += l*r; dF_1_dx += (1-y)*r -l*2*(1-x)*f1.f_x_00; dF_1_dy += -x* r +l*2*y* f1.f_y_11;
l = y* (1-x); r = ( f1.f_01 + x* x* f1.f_x_11 + (1-y)*(1-y)*f1.f_y_00 ); F_1 += l*r; dF_1_dx += -y* r +l*2*x* f1.f_x_11; dF_1_dy += (1-x)*r -l*2*(1-y)*f1.f_y_00;
l = y* x; r = ( f1.f_11 + (1-x)*(1-x)*f1.f_x_01 + (1-y)*(1-y)*f1.f_y_10 ); F_1 += l*r; dF_1_dx += y* r -l*2*(1-x)*f1.f_x_01; dF_1_dy += x* r -l*2*(1-y)*f1.f_y_10;
}
double result = (1-N_conj_ij)*F_0 + N_conj_ij*F_1;
*dFN_ij = (1-N_conj_ij)*dF_0_dx + N_conj_ij*dF_1_dx;
*dFN_ji = (1-N_conj_ij)*dF_0_dy + N_conj_ij*dF_1_dy;
*dFN_ij_conj = -F_0 + F_1;
return result;
}
/* ----------------------------------------------------------------------
read LCBOP potential file
------------------------------------------------------------------------- */
void PairLCBOP::read_file(char *filename)
{
int i,k,l;
char s[MAXLINE];
MPI_Comm_rank(world,&me);
// read file on proc 0
if (me == 0) {
FILE *fp = force->open_potential(filename);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open LCBOP potential file %s",filename);
error->one(FLERR,str);
}
// skip initial comment lines
while (1) {
fgets(s,MAXLINE,fp);
if (s[0] != '#') break;
}
// read parameters
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&r_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&r_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&gamma_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&A);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&B_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&B_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&alpha);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&beta_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&beta_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&d);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&C_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&C_4);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&C_6);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&L);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&kappa);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&R_0);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&R_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&r_0);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&r_1_LR);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&r_2_LR);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&v_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&v_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&eps_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&eps_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&lambda_1);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&lambda_2);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&eps);
fgets(s,MAXLINE,fp); sscanf(s,"%lg",&delta);
while (1) {
fgets(s,MAXLINE,fp);
if (s[0] != '#') break;
}
// F_conj spline
for (k = 0; k < 2; k++) { // 2 values of N_ij_conj
for (l = 0; l < 3; l++) { // 3 types of data: f, dfdx, dfdy
for (i = 0; i < 4; i++) { // 4x4 matrix
fgets(s,MAXLINE,fp);
sscanf(s,"%lg %lg %lg %lg",
&F_conj_data[i][0][k][l],
&F_conj_data[i][1][k][l],
&F_conj_data[i][2][k][l],
&F_conj_data[i][3][k][l]);
}
while (1) { fgets(s,MAXLINE,fp); if (s[0] != '#') break; }
}
}
// G spline
// x coordinates of mesh points
fgets(s,MAXLINE,fp);
sscanf( s,"%lg %lg %lg %lg %lg %lg",
&gX[0], &gX[1], &gX[2],
&gX[3], &gX[4], &gX[5] );
for (i = 0; i < 6; i++) { // for each power in polynomial
fgets(s,MAXLINE,fp);
sscanf( s,"%lg %lg %lg %lg %lg",
&gC[i][0], &gC[i][1], &gC[i][2],
&gC[i][3], &gC[i][4] );
}
fclose(fp);
}
// broadcast read-in and setup values
MPI_Bcast(&r_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&r_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&gamma_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&A ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&B_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&B_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&alpha ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&beta_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&beta_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&d ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&C_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&C_4 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&C_6 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&L ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&kappa ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&R_0 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&R_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&r_0 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&r_1_LR ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&r_2_LR ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&v_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&v_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&eps_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&eps_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&lambda_1 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&lambda_2 ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&eps ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&delta ,1,MPI_DOUBLE,0,world);
MPI_Bcast(&gX[0] ,6,MPI_DOUBLE,0,world);
MPI_Bcast(&gC[0][0] ,(6-1)*(5+1),MPI_DOUBLE,0,world);
MPI_Bcast(&F_conj_data[0],6*4*4,MPI_DOUBLE,0,world);
}
/* ----------------------------------------------------------------------
init coefficients for TF_conj
------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
#include <functional>
template< class function > void print_function( double x_0, double x_1, size_t n, function f, std::ostream &stream ) {
double dx = (x_1-x_0)/n;
for( double x=x_0; x<=x_1+0.0001; x+=dx ) {
double f_val, df;
f_val = f(x, &df);
stream << x << " " << f_val << " " << df << std::endl;
}
stream << std::endl;
}
void PairLCBOP::spline_init() {
for( size_t N_conj_ij=0; N_conj_ij<2; N_conj_ij++ ) // N_conj_ij
for( size_t N_ij=0; N_ij<4-1; N_ij++ )
for( size_t N_ji=0; N_ji<4-1; N_ji++ ) {
TF_conj_field &field = F_conj_field[N_ij][N_ji][N_conj_ij];
field.f_00 = F_conj_data[N_ij ][N_ji ][N_conj_ij][0];
field.f_01 = F_conj_data[N_ij ][N_ji+1][N_conj_ij][0];
field.f_10 = F_conj_data[N_ij+1][N_ji ][N_conj_ij][0];
field.f_11 = F_conj_data[N_ij+1][N_ji+1][N_conj_ij][0];
field.f_x_00 = F_conj_data[N_ij ][N_ji ][N_conj_ij][1] - field.f_10 + field.f_00;
field.f_x_01 = F_conj_data[N_ij ][N_ji+1][N_conj_ij][1] - field.f_11 + field.f_01;
field.f_x_10 = -(F_conj_data[N_ij+1][N_ji ][N_conj_ij][1] - field.f_10 + field.f_00);
field.f_x_11 = -(F_conj_data[N_ij+1][N_ji+1][N_conj_ij][1] - field.f_11 + field.f_01);
field.f_y_00 = F_conj_data[N_ij ][N_ji ][N_conj_ij][2] - field.f_01 + field.f_00;
field.f_y_01 = -(F_conj_data[N_ij ][N_ji+1][N_conj_ij][2] - field.f_01 + field.f_00);
field.f_y_10 = F_conj_data[N_ij+1][N_ji ][N_conj_ij][2] - field.f_11 + field.f_10;
field.f_y_11 = -(F_conj_data[N_ij+1][N_ji+1][N_conj_ij][2] - field.f_11 + field.f_10);
}
//some testing:
// std::ofstream file( "test.txt" );
// file << "gX:\n";
// file << gX[0] << " "
// << gX[1] << " "
// << gX[2] << " "
// << gX[3] << " "
// << gX[4] << " "
// << gX[5] << std::endl;
// file << "gC:\n";
// for( int i=0; i<6; i++ )
// file << gC[i][0] << " "
// << gC[i][1] << " "
// << gC[i][2] << " "
// << gC[i][3] << " "
// << gC[i][4] << std::endl;
// file << std::endl;
//
// file << "gamma_1 = " << gamma_1 << std::endl;
// file << "r_1 = " << r_1 << std::endl;
// file << "r_2 = " << r_2 << std::endl;
// file << "A = " << A << std::endl;
// file << "B_1 = " << B_1 << std::endl;
// file << "B_2 = " << B_2 << std::endl;
// file << "alpha = " << alpha << std::endl;
// file << "beta_1 = " << beta_1 << std::endl;
// file << "beta_2 = " << beta_2 << std::endl;
// file << "d = " << d << std::endl;
// file << "C_1 = " << C_1 << std::endl;
// file << "C_4 = " << C_4 << std::endl;
// file << "C_6 = " << C_6 << std::endl;
// file << "L = " << L << std::endl;
// file << "kappa = " << kappa << std::endl;
// file << "R_0 = " << R_0 << std::endl;
// file << "R_1 = " << R_1 << std::endl;
// file << "r_0 = " << r_0 << std::endl;
// file << "r_1_LR = " << r_1_LR << std::endl;
// file << "r_2_LR = " << r_2_LR << std::endl;
// file << "v_1 = " << v_1 << std::endl;
// file << "v_2 = " << v_2 << std::endl;
// file << "eps_1 = " << eps_1 << std::endl;
// file << "eps_2 = " << eps_2 << std::endl;
// file << "lambda_1 = " << lambda_1 << std::endl;
// file << "lambda_2 = " << lambda_2 << std::endl;
// file << "eps = " << eps << std::endl;
// file << "delta = " << delta << std::endl;
// file << "r_2_sq = " << r_2_sq << std::endl;
// file << std::endl;
//
//
// file << "gSpline:" << std::endl;
// double x_1 = 1, x_0 = -1;
// int n=1000;
// double dx = (x_1-x_0)/n;
// for( double x=x_0; x<=x_1+0.0001; x+=dx ) {
// double g, dg;
// g = gSpline(x, &dg);
// file << x << " " << g << " " << dg << std::endl;
// }
// file << std::endl;
//
// file << "hSpline:" << std::endl;
// double x_1 = 1, x_0 = -1;
// int n=1000;
// double dx = (x_1-x_0)/n;
// for( double x=x_0; x<=x_1+0.0001; x+=dx ) {
// double h, dh;
// h = hSpline(x, &dh);
// file << x << " " << h << " " << dh << std::endl;
// }
// file << std::endl;
//
//
// file << "f_c:" << std::endl;
// double x_1 = 4, x_0 = 0;
// int n=1000;
// double dx = (x_1-x_0)/n;
// for( double x=x_0; x<=x_1+0.0001; x+=dx ) {
// double f, df;
// f = f_c(x, r_1, r_2, &df);
// file << x << " " << f << " " << df << std::endl;
// }
// file << std::endl;
// file << "F_conj_data\n";
// for (int k = 0; k < 2; k++) { // 2 values of N_ij_conj
// for (int l = 0; l < 3; l++) { // 3 types of data: f, dfdx, dfdy
// for (int i = 0; i < 4; i++) { // 4x4 matrix
// file
// << F_conj_data[i][0][k][l] << " "
// << F_conj_data[i][1][k][l] << " "
// << F_conj_data[i][2][k][l] << " "
// << F_conj_data[i][3][k][l] << std::endl;
// }
// file << std::endl;
// }
// }
//
//
// file << "F_conj_0 ";
// double dummy;
// for( double y=0; y<=3.0+0.0001; y+=0.1 )
// file << y << " ";
// file << std::endl;
// for( double x=0; x<=3.0+0.0001; x+=0.1 ){
// file << x << " ";
// for( double y=0; y<=3.0+0.0001; y+=0.1 )
// file << F_conj( x, y, 0, &dummy, &dummy, &dummy ) << " ";
// file << std::endl;
// }
//
// file << "dF0_dx ";
// for( double y=0; y<=3.0+0.0001; y+=0.1 )
// file << y << " ";
// file << std::endl;
// for( double x=0; x<=3.0+0.0001; x+=0.1 ){
// file << x << " ";
// for( double y=0; y<=3.0+0.0001; y+=0.1 ) {
// double dF_dx;
// F_conj( x, y, 0, &dF_dx, &dummy, &dummy );
// file << dF_dx << " ";
// }
// file << std::endl;
// }
//
//
//
// file << "F_conj_1 ";
// for( double y=0; y<=3.0+0.0001; y+=0.1 )
// file << y << " ";
// file << std::endl;
// for( double x=0; x<=3.0+0.0001; x+=0.1 ){
// file << x << " ";
// for( double y=0; y<=3.0+0.0001; y+=0.1 )
// file << F_conj( x, y, 0, &dummy, &dummy, &dummy ) << " ";
// file << std::endl;
// }
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairLCBOP::memory_usage()
{
double bytes = 0.0;
bytes += maxlocal * sizeof(int);
bytes += maxlocal * sizeof(int *);
for (int i = 0; i < comm->nthreads; i++)
bytes += ipage[i].size();
bytes += 3*maxlocal * sizeof(double);
return bytes;
}
Event Timeline
Log In to Comment