diff --git a/src/MANYBODY/pair_lcbop.cpp b/src/MANYBODY/pair_lcbop.cpp index 97c1bfef0b..691f1f0b11 100644 --- a/src/MANYBODY/pair_lcbop.cpp +++ b/src/MANYBODY/pair_lcbop.cpp @@ -1,1295 +1,1295 @@ -/* ---------------------------------------------------------------------- - 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 "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; - ghostneigh = 1; - - maxlocal = 0; - SR_numneigh = NULL; - SR_firstneigh = NULL; - maxpage = 0; - pages = NULL; - N = NULL; - M = NULL; -} - -/* ---------------------------------------------------------------------- - check if allocated, since class can be destructed when incomplete -------------------------------------------------------------------------- */ - -PairLCBOP::~PairLCBOP() { - memory->destroy(SR_numneigh); - memory->sfree(SR_firstneigh); - for (int i = 0; i < maxpage; i++) memory->destroy(pages[i]); - memory->sfree(pages); - 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); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; - neighbor->requests[irequest]->ghost = 1; - - // local SR neighbor list memory - - pgsize = neighbor->pgsize; - oneatom = neighbor->oneatom; - if (maxpage == 0) add_pages(); -} - -/* ---------------------------------------------------------------------- - 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, // number of atoms(both local and ghost) neighbors are stored for - jnum; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; - int *ilist,*jlist,*numneigh,**firstneigh; - int *neighptr; - - double **x = atom->x; - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - - 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 - - int npage = 0; - int npnt = 0; // position in current page - - for (ii = 0; ii < allnum; ii++) { - i = ilist[ii]; - - if (pgsize - npnt < oneatom) { // ensure at least oneatom space free at current page - npnt = 0; - npage++; - if (npage == maxpage) add_pages(); - } - neighptr = &pages[npage][npnt]; - n = 0; - - 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; - npnt += n; - if( npnt >= pgsize ) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); - } - - - // 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 Nij = N[i]-f_c_ij; - // F(xij) = 1-f_c_LR(Nji, 2,3,&dummy) - M[i] += f_c_ij * ( 1-f_c_LR(Nij, 2,3,&dS) ); - } - } - } -} - -/* ---------------------------------------------------------------------- - Short range forces and energy -------------------------------------------------------------------------- */ - -void PairLCBOP::FSR(int eflag, int vflag) -{ - int i,j,jj,ii,inum,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 d_f_c_ij,del[3]; - int *ilist,*SR_neighs; - - double **x = atom->x; - double **f = atom->f; - int *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,m,ii,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; - int *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 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 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; - double f_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 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; -} - -/* ---------------------------------------------------------------------- - add pages to SR neighbor list -------------------------------------------------------------------------- */ - -void PairLCBOP::add_pages(int howmany) -{ - int toppage = maxpage; - maxpage += howmany*PGDELTA; - - pages = (int **) - memory->srealloc(pages,maxpage*sizeof(int *),"LCBOP:pages"); - for (int i = toppage; i < maxpage; i++) - memory->create(pages[i],pgsize,"LCBOP:pages[i]"); -} - -/* ---------------------------------------------------------------------- - 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 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( floor( N_ij ) ), 2 ); // 2 is the highest number of field - size_t N_ji_int = MIN( static_cast( 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,j,k,l,limit; - char s[MAXLINE]; - - MPI_Comm_rank(world,&me); - - // read file on proc 0 - - if (me == 0) { - FILE *fp = fopen(filename,"r"); - 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 -#include -#include -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 *); - bytes += maxpage * neighbor->pgsize * sizeof(int); - bytes += 3 * maxlocal * sizeof(double); - return bytes; -} +/* ---------------------------------------------------------------------- + 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 "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; + ghostneigh = 1; + + maxlocal = 0; + SR_numneigh = NULL; + SR_firstneigh = NULL; + maxpage = 0; + pages = NULL; + N = NULL; + M = NULL; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairLCBOP::~PairLCBOP() { + memory->destroy(SR_numneigh); + memory->sfree(SR_firstneigh); + for (int i = 0; i < maxpage; i++) memory->destroy(pages[i]); + memory->sfree(pages); + 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); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->ghost = 1; + + // local SR neighbor list memory + + pgsize = neighbor->pgsize; + oneatom = neighbor->oneatom; + if (maxpage == 0) add_pages(); +} + +/* ---------------------------------------------------------------------- + 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, // number of atoms(both local and ghost) neighbors are stored for + jnum; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; + int *ilist,*jlist,*numneigh,**firstneigh; + int *neighptr; + + double **x = atom->x; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + 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 + + int npage = 0; + int npnt = 0; // position in current page + + for (ii = 0; ii < allnum; ii++) { + i = ilist[ii]; + + if (pgsize - npnt < oneatom) { // ensure at least oneatom space free at current page + npnt = 0; + npage++; + if (npage == maxpage) add_pages(); + } + neighptr = &pages[npage][npnt]; + n = 0; + + 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; + npnt += n; + if( npnt >= pgsize ) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + + + // 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,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 d_f_c_ij,del[3]; + int *ilist,*SR_neighs; + + double **x = atom->x; + double **f = atom->f; + int *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,m,ii,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; + int *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 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 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; + double f_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 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; +} + +/* ---------------------------------------------------------------------- + add pages to SR neighbor list +------------------------------------------------------------------------- */ + +void PairLCBOP::add_pages(int howmany) +{ + int toppage = maxpage; + maxpage += howmany*PGDELTA; + + pages = (int **) + memory->srealloc(pages,maxpage*sizeof(int *),"LCBOP:pages"); + for (int i = toppage; i < maxpage; i++) + memory->create(pages[i],pgsize,"LCBOP:pages[i]"); +} + +/* ---------------------------------------------------------------------- + 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 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( floor( N_ij ) ), 2 ); // 2 is the highest number of field + size_t N_ji_int = MIN( static_cast( 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,j,k,l,limit; + char s[MAXLINE]; + + MPI_Comm_rank(world,&me); + + // read file on proc 0 + + if (me == 0) { + FILE *fp = fopen(filename,"r"); + 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 +#include +#include +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 *); + bytes += maxpage * neighbor->pgsize * sizeof(int); + bytes += 3 * maxlocal * sizeof(double); + return bytes; +}