From b8c97b48aff13c2af4ad8539b0df1d4bd2b0b3fc Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 11 Feb 2009 22:26:01 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2568 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/REAX/Install.csh | 32 +-- src/REAX/pair_reax.cpp | 497 +++++++++++++++-------------------- src/REAX/pair_reax.h | 15 +- src/REAX/pair_reax_fortran.h | 201 ++++++++++++++ src/REAX/reax_cbkabo.h | 16 -- src/REAX/reax_cbkbo.h | 4 - src/REAX/reax_cbkc.h | 6 - src/REAX/reax_cbkch.h | 10 - src/REAX/reax_cbkd.h | 12 - src/REAX/reax_cbkia.h | 6 - src/REAX/reax_cbklonpar.h | 5 - src/REAX/reax_cbknubon2.h | 6 - src/REAX/reax_cbkpairs.h | 15 -- src/REAX/reax_cbkqa.h | 5 - src/REAX/reax_energies.h | 25 -- src/REAX/reax_fortran.h | 31 --- src/REAX/reax_functions.h | 37 --- src/REAX/reax_params.h | 14 - src/REAX/reax_small.h | 15 -- 19 files changed, 420 insertions(+), 532 deletions(-) create mode 100644 src/REAX/pair_reax_fortran.h delete mode 100644 src/REAX/reax_cbkabo.h delete mode 100644 src/REAX/reax_cbkbo.h delete mode 100644 src/REAX/reax_cbkc.h delete mode 100644 src/REAX/reax_cbkch.h delete mode 100644 src/REAX/reax_cbkd.h delete mode 100644 src/REAX/reax_cbkia.h delete mode 100644 src/REAX/reax_cbklonpar.h delete mode 100644 src/REAX/reax_cbknubon2.h delete mode 100644 src/REAX/reax_cbkpairs.h delete mode 100644 src/REAX/reax_cbkqa.h delete mode 100644 src/REAX/reax_energies.h delete mode 100644 src/REAX/reax_fortran.h delete mode 100644 src/REAX/reax_functions.h delete mode 100644 src/REAX/reax_params.h delete mode 100644 src/REAX/reax_small.h diff --git a/src/REAX/Install.csh b/src/REAX/Install.csh index 3e587f39fb..ef6e503691 100644 --- a/src/REAX/Install.csh +++ b/src/REAX/Install.csh @@ -7,21 +7,7 @@ if ($1 == 1) then cp pair_reax.cpp .. cp pair_reax.h .. - cp reax_cbkabo.h .. - cp reax_cbkbo.h .. - cp reax_cbkc.h .. - cp reax_cbkch.h .. - cp reax_cbkd.h .. - cp reax_cbkia.h .. - cp reax_cbklonpar.h .. - cp reax_cbknubon2.h .. - cp reax_cbkpairs.h .. - cp reax_cbkqa.h .. - cp reax_energies.h .. - cp reax_fortran.h .. - cp reax_functions.h .. - cp reax_params.h .. - cp reax_small.h .. + cp pair_reax_fortran.h .. else if ($1 == 0) then @@ -31,20 +17,6 @@ else if ($1 == 0) then rm ../pair_reax.cpp rm ../pair_reax.h - rm ../reax_cbkabo.h - rm ../reax_cbkbo.h - rm ../reax_cbkc.h - rm ../reax_cbkch.h - rm ../reax_cbkd.h - rm ../reax_cbkia.h - rm ../reax_cbklonpar.h - rm ../reax_cbknubon2.h - rm ../reax_cbkpairs.h - rm ../reax_cbkqa.h - rm ../reax_energies.h - rm ../reax_fortran.h - rm ../reax_functions.h - rm ../reax_params.h - rm ../reax_small.h + rm ../pair_reax_fortran.h endif diff --git a/src/REAX/pair_reax.cpp b/src/REAX/pair_reax.cpp index 679635f3f6..b10d570b72 100644 --- a/src/REAX/pair_reax.cpp +++ b/src/REAX/pair_reax.cpp @@ -12,8 +12,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Hansohl Cho (MIT, hansohl@mit.edu) - Aidan Thompson (Sandia, athomps@sandia.gov) + Contributing authors: Aidan Thompson (Sandia, athomps@sandia.gov) + Hansohl Cho (MIT, hansohl@mit.edu) LAMMPS implementation of the Reactive Force Field (ReaxFF) is based on Aidan Thompson's GRASP code (General Reactive Atomistic Simulation Program) @@ -29,30 +29,17 @@ #include "stdio.h" #include "stdlib.h" #include "string.h" - #include "pair_reax.h" +#include "pair_reax_fortran.h" #include "atom.h" #include "force.h" #include "comm.h" -#include "memory.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "memory.h" #include "error.h" -#include "reax_fortran.h" -#include "reax_params.h" -#include "reax_cbkc.h" -#include "reax_cbkd.h" -#include "reax_cbkch.h" -#include "reax_cbkabo.h" -#include "reax_cbkia.h" -#include "reax_cbkpairs.h" -#include "reax_energies.h" -#include "reax_small.h" -#include "reax_functions.h" - using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) @@ -67,14 +54,20 @@ PairREAX::PairREAX(LAMMPS *lmp) : Pair(lmp) one_coeff = 1; Lmidpoint = false; - cutmax=0.0; - hbcut=6.0; - iprune=4; - ihb=1; - + cutmax = 0.0; + hbcut = 6.0; + iprune = 4; + ihb = 1; chpot = 0; - comm_forward = 2; + nmax = 0; + rcg = NULL; + wcg = NULL; + pcg = NULL; + poldcg = NULL; + qcg = NULL; + + comm_forward = 1; comm_reverse = 1; precision = 1.0e-6; @@ -90,44 +83,59 @@ PairREAX::~PairREAX() if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); + + for (int i = 1; i <= atom->ntypes; i++) + delete [] param_list[i].params; + delete [] param_list; } + + memory->sfree(rcg); + memory->sfree(wcg); + memory->sfree(pcg); + memory->sfree(poldcg); + memory->sfree(qcg); } /* ---------------------------------------------------------------------- */ void PairREAX::compute(int eflag, int vflag) { + int i,j; double evdwl,ecoul; - double reax_energy_pieces[13]; - double energy_charge_equilibration; - evdwl = 0.0; - ecoul = 0.0; - + evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - int newton_pair = force->newton_pair; - - if (eflag) - for (int k = 0; k < 14; k++) - reax_energy_pieces[k]=0; - - if (vflag) FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 1; + if (vflag_global) FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 1; else FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 0; - if (evflag) { - FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 1; - FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 1; - } else FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 0; + if (vflag_atom) FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 1; + else FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 0; + + // reallocate CG arrays if necessary + + if (atom->nmax > nmax) { + memory->sfree(rcg); + memory->sfree(wcg); + memory->sfree(pcg); + memory->sfree(poldcg); + memory->sfree(qcg); + nmax = atom->nmax; + int n = nmax+1; + rcg = (double *) memory->smalloc(n*sizeof(double),"reax:rcg"); + wcg = (double *) memory->smalloc(n*sizeof(double),"reax:wcg"); + pcg = (double *) memory->smalloc(n*sizeof(double),"reax:pcg"); + poldcg = (double *) memory->smalloc(n*sizeof(double),"reax:poldcg"); + qcg = (double *) memory->smalloc(n*sizeof(double),"reax:qcg"); + } // calculate the atomic charge distribution - + compute_charge(energy_charge_equilibration); - // transfer LAMMPS positions and neighbor lists - // to REAX Fortran positions and Verlet lists + // transfer LAMMPS positions and neighbor lists to REAX write_reax_positions(); write_reax_vlist(); @@ -142,75 +150,49 @@ void PairREAX::compute(int eflag, int vflag) // communicate local atomic bond order to ghost atomic bond order + packflag = 0; comm->comm_pair(this); FORTRAN(molec, MOLEC)(); FORTRAN(encalc, ENCALC)(); - - int node = comm -> me; - FORTRAN(mdsav, MDSAV)(&node); + FORTRAN(mdsav, MDSAV)(&comm->me); // read forces from ReaxFF Fortran read_reax_forces(); - // read energy from ReaxFF Fortran + // extract global and per-atom energy from ReaxFF Fortran - if (eflag) { - reax_energy_pieces[0] = FORTRAN(cbkenergies, CBKENERGIES).eb; - reax_energy_pieces[1] = FORTRAN(cbkenergies, CBKENERGIES).ea; - reax_energy_pieces[2] = FORTRAN(cbkenergies, CBKENERGIES).elp; - reax_energy_pieces[3] = FORTRAN(cbkenergies, CBKENERGIES).emol; - reax_energy_pieces[4] = FORTRAN(cbkenergies, CBKENERGIES).ev; - reax_energy_pieces[5] = FORTRAN(cbkenergies, CBKENERGIES).epen; - reax_energy_pieces[6] = FORTRAN(cbkenergies, CBKENERGIES).ecoa; - reax_energy_pieces[7] = FORTRAN(cbkenergies, CBKENERGIES).ehb; - reax_energy_pieces[8] = FORTRAN(cbkenergies, CBKENERGIES).et; - reax_energy_pieces[9] = FORTRAN(cbkenergies, CBKENERGIES).eco; - reax_energy_pieces[10] = FORTRAN(cbkenergies, CBKENERGIES).ew; - reax_energy_pieces[11] = FORTRAN(cbkenergies, CBKENERGIES).ep; - reax_energy_pieces[12] = FORTRAN(cbkenergies, CBKENERGIES).efi; - - for (int k = 0; k < 13; k++) - evdwl += reax_energy_pieces[k]; + if (eflag_global) { + evdwl += FORTRAN(cbkenergies, CBKENERGIES).eb; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).ea; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).elp; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).emol; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).ev; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).epen; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).ecoa; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).ehb; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).et; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).eco; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).ew; + evdwl += FORTRAN(cbkenergies, CBKENERGIES).efi; - // LAMMPS eVDWL energy does not include the Coulomb energy in ReaxFF - - evdwl = evdwl - reax_energy_pieces[11]; - - // LAMMPS eCOUL energy includes the Coulomb energy and - // charge equilibation energy - - ecoul = reax_energy_pieces[11] + energy_charge_equilibration; + ecoul += FORTRAN(cbkenergies, CBKENERGIES).ep; + ecoul += energy_charge_equilibration; eng_vdwl += evdwl; eng_coul += ecoul; } - int nval, ntor, nhb, nbonall, nbon, na, na_local, nvpair; - int reaxsize[8]; + if (eflag_atom) { + int ntotal = atom->nlocal + atom->nghost; + for (i = 0; i < ntotal; i++) + eatom[i] = FORTRAN(cbkd,CBKD).estrain[i]; + } - na_local = FORTRAN(rsmall, RSMALL).na_local; - na = FORTRAN(rsmall, RSMALL).na; - FORTRAN(getnbonall, GETNBONALL)(&nbonall); - nbon = FORTRAN(rsmall, RSMALL).nbon; - FORTRAN(getnval, GETNVAL)(&nval); - FORTRAN(getntor, GETNTOR)(&ntor); - FORTRAN(getnhb, GETNHB)(&nhb); - nvpair = FORTRAN(cbkpairs, CBKPAIRS).nvpair; + // extract global and per-atom virial from ReaxFF Fortran - reaxsize[0] = na_local; - reaxsize[1] = na; - reaxsize[2] = nbonall; - reaxsize[3] = nbon; - reaxsize[4] = nval; - reaxsize[5] = ntor; - reaxsize[6] = nhb; - reaxsize[7] = nvpair; - - // call ev_tally to update the global energy and virial - - if (vflag) { + if (vflag_global) { virial[0] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[0]; virial[1] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[1]; virial[2] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[2]; @@ -219,7 +201,19 @@ void PairREAX::compute(int eflag, int vflag) virial[5] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[5]; } - if (vflag_atom) read_reax_atom_virial(); + if (vflag_atom) { + int ntotal = atom->nlocal + atom->nghost; + j = 0; + for (i = 0; i < ntotal; i++) { + vatom[i][0] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+0]; + vatom[i][1] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+1]; + vatom[i][2] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+2]; + vatom[i][3] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+3]; + vatom[i][4] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+4]; + vatom[i][5] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+5]; + j += 6; + } + } } /* ---------------------------------------------------------------------- */ @@ -343,22 +337,18 @@ void PairREAX::write_reax_vlist() FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 0; - // this is for Lmidpoint = false - - if (i < nlocal) { - if (j < nlocal) + if (j < nlocal) + FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; + else if (itag < jtag) + FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; + else if (itag == jtag) { + if (delz > SMALL) FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (itag < jtag) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (itag == jtag) { - if (delz > SMALL) + else if (fabs(delz) < SMALL) { + if (dely > SMALL) + FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; + else if (fabs(dely) < SMALL && delx > SMALL) FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (fabs(delz) < SMALL) { - if (dely > SMALL) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (fabs(dely) < SMALL && delx > SMALL) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - } } } nvpair++; @@ -406,25 +396,6 @@ void PairREAX::write_reax_vlist() } FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 0; - - // this is for Lmidpoint = false - - if (i < nlocal) { - if (j < nlocal) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (itag < jtag) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (itag == jtag) { - if (delz > SMALL) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (fabs(delz) < SMALL) { - if (dely > SMALL) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - else if (fabs(dely) < SMALL && delx > SMALL) - FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1; - } - } - } nvpair++; } } @@ -461,34 +432,6 @@ void PairREAX::read_reax_forces() /* ---------------------------------------------------------------------- */ -void PairREAX::read_reax_atom_virial() -{ - double vatomtmp[6]; - - int ntotal = atom->nlocal + atom->nghost; - - int j = 0; - for (int i =0; i < ntotal; i++) { - vatomtmp[0] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+0]; - vatomtmp[1] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+1]; - vatomtmp[2] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+2]; - vatomtmp[3] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+3]; - vatomtmp[4] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+4]; - vatomtmp[5] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+5]; - - vatom[i][0] = vatomtmp[0]; - vatom[i][1] = vatomtmp[1]; - vatom[i][2] = vatomtmp[2]; - vatom[i][3] = vatomtmp[3]; - vatom[i][4] = vatomtmp[4]; - vatom[i][5] = vatomtmp[5]; - - j += 6; - } -} - -/* ---------------------------------------------------------------------- */ - void PairREAX::allocate() { allocated = 1; @@ -496,6 +439,10 @@ void PairREAX::allocate() setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + param_list = new ff_params[n+1]; + for (int i = 1; i <= n; i++) + param_list[i].params = new double[5]; } /* ---------------------------------------------------------------------- @@ -549,13 +496,13 @@ void PairREAX::init_style() { // REAX requires pair newton off + if (atom->tag_enable == 0) + error->all("Pair style reax requires atom IDs"); if (force->newton_pair) - error->all("Pair style reax requires newton pair off"); + error->all("Pair style reax requires newton pair on"); int irequest = neighbor->request(this); - int node = comm->me; - FORTRAN(readc, READC)(); FORTRAN(init, INIT)(); FORTRAN(ffinpt, FFINPT)(); @@ -565,7 +512,7 @@ void PairREAX::init_style() int ngeofor_tmp = -1; FORTRAN(setngeofor, SETNGEOFOR)(&ngeofor_tmp); - if (node == 0) FORTRAN(readgeo, READGEO)(); + if (comm->me == 0) FORTRAN(readgeo, READGEO)(); // initial setup for cutoff radius of VLIST and BLIST in ReaxFF @@ -580,23 +527,18 @@ void PairREAX::init_style() // parameters for charge equilibration from ReaxFF input, fort.4 FORTRAN(getswa, GETSWA)(&swa); - ff_params ff_param_tmp; double chi, eta, gamma; - int nparams = 5; - param_list = new ff_params[atom->ntypes+1]; for (int itype = 1; itype <= atom->ntypes; itype++) { chi = FORTRAN(cbkchb, CBKCHB).chi[itype-1]; eta = FORTRAN(cbkchb, CBKCHB).eta[itype-1]; gamma = FORTRAN(cbkchb, CBKCHB).gam[itype-1]; - ff_param_tmp.np = nparams; - ff_param_tmp.rcutsq = cutmax; - ff_param_tmp.params = new double[nparams]; - ff_param_tmp.params[0] = chi; - ff_param_tmp.params[1] = eta; - ff_param_tmp.params[2] = gamma; - ff_param_tmp.params[3] = swa; - ff_param_tmp.params[4] = swb; - param_list[itype] = ff_param_tmp; + param_list[itype].np = 5; + param_list[itype].rcutsq = cutmax; + param_list[itype].params[0] = chi; + param_list[itype].params[1] = eta; + param_list[itype].params[2] = gamma; + param_list[itype].params[3] = swa; + param_list[itype].params[4] = swb; } taper_setup(); @@ -618,12 +560,21 @@ int PairREAX::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) int i,j,m; m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = FORTRAN(cbkabo, CBKABO).abo[j]; - buf[m++] = w[j]; + + if (packflag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = FORTRAN(cbkabo, CBKABO).abo[j]; + } + + } else { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = wcg[j]; + } } - return 2; + + return 1; } /* ---------------------------------------------------------------------- */ @@ -634,9 +585,14 @@ void PairREAX::unpack_comm(int n, int first, double *buf) m = 0; last = first + n; - for (i = first; i < last; i++) { - FORTRAN(cbkabo, CBKABO).abo[i] = buf[m++]; - w[i] = buf[m++]; + + if (packflag == 0) { + for (i = first; i < last; i++) + FORTRAN(cbkabo, CBKABO).abo[i] = buf[m++]; + + } else { + for (i = first; i < last; i++) + wcg[i] = buf[m++]; } } @@ -649,7 +605,7 @@ int PairREAX::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; for (i = first; i < last; i++) - buf[m++] = w[i]; + buf[m++] = wcg[i]; return 1; } @@ -663,7 +619,7 @@ void PairREAX::unpack_reverse_comm(int n, int *list, double *buf) m = 0; for (i = 0; i < n; i++) { j = list[i]; - w[j] += buf[m++]; + wcg[j] += buf[m++]; } } @@ -721,14 +677,10 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) { double xitmp, yitmp, zitmp; double xjtmp, yjtmp, zjtmp; - int itype, jtype, itag, jtag; - int ii, jj, i, j; double delr2, delr_norm, gamt, hulp1, hulp2; double delx, dely, delz; - int node, nprocs; - double qsum, qi; double *ch; double *elcvec; @@ -738,6 +690,10 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) int maxmatentries, nmatentries; double sw; double rtmp[3]; + int inum, jnum; + int *ilist; + int *jlist; + int *numneigh, **firstneigh; double **x = atom->x; double *q = atom->q; @@ -747,19 +703,13 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) int nlocal = atom->nlocal; int nghost = atom->nghost; - int inum, jnum; - int *ilist; - int *jlist; - int *numneigh, **firstneigh; - - node = comm -> me; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; int numneigh_total = 0; - for (ii =0; ii < inum; ii++) + for (ii = 0; ii < inum; ii++) numneigh_total += numneigh[ilist[ii]]; arow_ptr = new int[nlocal+nghost+1]; @@ -772,7 +722,7 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) // build linear system - for (ii =0; ii= nlocal) { - if (neighbor->style==0) - { - if (itag > jtag) - { - if ((itag+jtag) % 2 == 0) continue; - } - else if (itag < jtag) - { - if ((itag+jtag) % 2 == 1) continue; - } - else - // if (itag == jtag) - { - if (zjtmp < zitmp) continue; - else if (zjtmp == zitmp && yjtmp < yitmp) continue; - else if (zjtmp == zitmp && yjtmp == yitmp && xjtmp < xitmp) - continue; - } - } - else if (neighbor->style==1) - { - if (zjtmp < zitmp) continue; - else if (zjtmp == zitmp && yjtmp < yitmp) continue; - else if (zjtmp == zitmp && yjtmp == yitmp && xjtmp < xitmp) - continue; - } - else - { - error->all("ReaxFF does not support multi style neiborlist"); - } - } - - - // rcutvsq = cutmax*cutmax, in ReaxFF - if (delr2 <= rcutvsq) - { - gamt = sqrt(param_list[itype].params[2]* - param_list[jtype].params[2]); - delr_norm = sqrt(delr2); - sw = taper_E(delr_norm, delr2); - hulp1=(delr_norm*delr2+(1.0/(gamt*gamt*gamt))); - hulp2=sw*14.40/cbrt(hulp1); - aval[nmatentries] = hulp2; - acol_ind[nmatentries] = j; - nmatentries++; + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (zjtmp < zitmp) continue; + if (zjtmp == zitmp && yjtmp < yitmp) continue; + if (zjtmp == zitmp && yjtmp == yitmp && xjtmp < xitmp) continue; } + } + + // rcutvsq = cutmax*cutmax, in ReaxFF + + if (delr2 <= rcutvsq) { + gamt = sqrt(param_list[itype].params[2]*param_list[jtype].params[2]); + delr_norm = sqrt(delr2); + sw = taper_E(delr_norm, delr2); + hulp1=(delr_norm*delr2+(1.0/(gamt*gamt*gamt))); + hulp2=sw*14.40/cbrt(hulp1); + aval[nmatentries] = hulp2; + acol_ind[nmatentries] = j; + nmatentries++; + } } } - // In this case, we don't use Midpoint method - // So, we don't need to consider ghost-ghost interactions - // But, need to fill the arow_ptr[] arrays for the ghost atoms + // in this case, we don't use Midpoint method + // so, we don't need to consider ghost-ghost interactions + // but, need to fill the arow_ptr[] arrays for the ghost atoms for (i = nlocal; i < nlocal+nghost; i++) arow_ptr[i] = nmatentries; @@ -877,14 +803,14 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) for (i = nlocal; i < nlocal+nghost; i++) elcvec[i] = 0.0; - // Assign current charges to charge vector + // assign current charges to charge vector qsum = 0.0; - for (ii =0; iireverse_comm_pair(this); comm->comm_pair(this); @@ -1044,45 +962,30 @@ void PairREAX::cg_solve(const int & nlocal, const int & nghost, for (int i=0; i Fortran calling syntax +// CONS(a,b) should return ab, the concatenation of its arguments + +#if __STDC__ +#define CONS(a,b) a##b +#else +#define CONS(a,b) a/**/b +#endif + +#ifdef _ARDENT +#define FORTRAN(lcname,ucname) ucname +#endif + +#ifdef _IBM +#define FORTRAN(lcname,ucname) lcname +#endif + +#ifdef _HP +#define FORTRAN(lcname,ucname) lcname +#endif + +#ifdef _F2C_LINUX +#define FORTRAN(lcname,ucname) CONS(lcname,__) +#endif + +#ifndef FORTRAN +#define FORTRAN(lcname,ucname) CONS(lcname,_) +#endif + +// hard-wired array sizes set in Fortran library + +#include "reax_defs.h" + +class ReaxParams { + public: + enum {nneighmax=NNEIGHMAXDEF, + nat=NATDEF, + nattot=NATTOTDEF, + nsort=NSORTDEF, + mbond=MBONDDEF, + nbomax=NBOMAXDEF, + }; +}; + +// data structures corresponding to values in Fortran library + +extern struct { + double abo[ReaxParams::nat]; +} FORTRAN(cbkabo,CBKABO); + +extern struct { + double bo[ReaxParams::nbomax]; +} FORTRAN(cbkbo,CBKBO); + +extern struct { + double c[3*ReaxParams::nat]; double cglobal[3*ReaxParams::nattot]; + int itag[ReaxParams::nat]; +} FORTRAN(cbkc,CBKC); + +extern struct { +double ch[ReaxParams::nat]; +} FORTRAN(cbkch,CBKCH); + +extern struct { + double chi[ReaxParams::nsort]; + double eta[ReaxParams::nsort]; + double gam[ReaxParams::nsort]; +} FORTRAN(cbkchb,CBKCHB); + +extern struct { + double d[3*ReaxParams::nat]; double estrain[ReaxParams::nat]; +} FORTRAN(cbkd,CBKD); + +extern struct { + double atomvirial[6*ReaxParams::nat]; + double virial[6]; + int Lvirial; + int Latomvirial; +} FORTRAN(cbkvirial,CBKVIRIAL); + +extern struct { + int ia[ReaxParams::nat*(ReaxParams::mbond+3)]; + int iag[ReaxParams::nat*(ReaxParams::mbond+3)]; +} FORTRAN(cbkia,CBKIA); + +extern struct { + double vlp[ReaxParams::nat]; + double dvlpdsbo[ReaxParams::nat]; +} FORTRAN(cbklonpar,CBKLONPAR); + +extern struct { + int nubon1[ReaxParams::nat*(ReaxParams::mbond)]; + int nubon2[ReaxParams::nat*(ReaxParams::mbond)]; +} FORTRAN(cbknubon2,CBKNUBON2); + +extern struct { + int nvl1[ReaxParams::nneighmax * ReaxParams::nat]; + int nvl2[ReaxParams::nneighmax * ReaxParams::nat]; + int nvpair; + int nvlself; +} FORTRAN(cbkpairs,CBKPAIRS); + +extern struct { + int nvlbo[ReaxParams::nneighmax * ReaxParams::nat]; +} FORTRAN(cbknvlbo,CBKNVLBO); + +extern struct { + int nvlown[ReaxParams::nneighmax * ReaxParams::nat]; +} FORTRAN(cbknvlown,CBKNVLOWN); + +extern struct { + char qa[20*ReaxParams::nattot+10]; +} FORTRAN(cbkqa,CBKQA); + +extern struct { + double eb; + double eoop; + double epen; + double estrc; + double deda[3]; + double pressu; + double efi; + double elp; + double emol; + double ea; + double eres; + double et; + double eradbo; + double ev; + double eco; + double ecoa; + double ehb; + double sw; + double ew; + double ep; + double ekin; +} FORTRAN(cbkenergies,CBKENERGIES); + +extern struct { + double tset; + double dseed; + double tempmd; + double ts2; + double ts22; + int nmolo; + int nmolo5; + int nbon; + int na; + int namov; + int na_local; +} FORTRAN(rsmall,RSMALL); + +// external routines provided by Fortran library + +extern "C" void FORTRAN(readc,READC)(); +extern "C" void FORTRAN(init,INIT)(); +extern "C" void FORTRAN(ffinpt,FFINPT)(); +extern "C" void FORTRAN(tap7th,TAP7TH)(); +extern "C" void FORTRAN(taper,TAPER)(double*,double*); +extern "C" void FORTRAN(readgeo,READGEO)(); +extern "C" void FORTRAN(srtatom,SRTATOM)(); +extern "C" void FORTRAN(vlist,VLIST) (); +extern "C" void FORTRAN(srtbon1,SRTBON1)(int*,int*,double*); +extern "C" void FORTRAN(molec,MOLEC)(); +extern "C" void FORTRAN(encalc,ENCALC)(); +extern "C" void FORTRAN(getswb,GETSWB)(double*); +extern "C" void FORTRAN(getswa,GETSWA)(double*); +extern "C" void FORTRAN(getvrange,GET_VRANGE)(double*); +extern "C" void FORTRAN(getnvlist,GET_NVLIST)(int*); +extern "C" void FORTRAN(getvlbora,GETVLBORA)(double*); +extern "C" void FORTRAN(cgsolve,CGSOLVE) + (int*,double*,int*,double*,double*,int*); +extern "C" void FORTRAN(getnval,GETNVAL)(int*); +extern "C" void FORTRAN(getntor,GETNTOR)(int*); +extern "C" void FORTRAN(getnhb,GETNHB)(int*); +extern "C" void FORTRAN(getnbonall,GETNBONALL)(int*); +extern "C" void FORTRAN(getnneighmax,GETNNEIGHMAX)(int*); +extern "C" void FORTRAN(getnat,GETNAT)(int*); +extern "C" void FORTRAN(getnattot,GETNATTOT)(int*); +extern "C" void FORTRAN(getnsort,GETNSORT)(int*); +extern "C" void FORTRAN(getmbond,GETMBOND)(int*); +extern "C" void FORTRAN(getnso,GETNSO)(int*); +extern "C" void FORTRAN(setngeofor,SETNGEOFOR)(int*); +extern "C" void FORTRAN(mdsav,MDSAV)(int*); +extern "C" void FORTRAN(getnsbmax,GETNSBMAX)(int*); +extern "C" void FORTRAN(getnsbma2,GETNSBMA2)(int*); +extern "C" void FORTRAN(getcutof3,GETCUTOF3)(double*); + diff --git a/src/REAX/reax_cbkabo.h b/src/REAX/reax_cbkabo.h deleted file mode 100644 index ab725b571a..0000000000 --- a/src/REAX/reax_cbkabo.h +++ /dev/null @@ -1,16 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -extern struct { - double abo[ReaxParams::nat]; -} FORTRAN(cbkabo,CBKABO); diff --git a/src/REAX/reax_cbkbo.h b/src/REAX/reax_cbkbo.h deleted file mode 100644 index b276a9084a..0000000000 --- a/src/REAX/reax_cbkbo.h +++ /dev/null @@ -1,4 +0,0 @@ -/////:EOH~ -extern struct { - double bo[ReaxParams::nbomax]; -} FORTRAN(cbkbo,CBKBO); diff --git a/src/REAX/reax_cbkc.h b/src/REAX/reax_cbkc.h deleted file mode 100644 index cbb275a459..0000000000 --- a/src/REAX/reax_cbkc.h +++ /dev/null @@ -1,6 +0,0 @@ -/////:EOH~ -extern struct { - double c[3*ReaxParams::nat]; double cglobal[3*ReaxParams::nattot]; - int itag[ReaxParams::nat]; -} FORTRAN(cbkc,CBKC); - diff --git a/src/REAX/reax_cbkch.h b/src/REAX/reax_cbkch.h deleted file mode 100644 index f25b9fb46b..0000000000 --- a/src/REAX/reax_cbkch.h +++ /dev/null @@ -1,10 +0,0 @@ -/////:EOH~ -extern struct { -double ch[ReaxParams::nat]; -} FORTRAN(cbkch,CBKCH); - -extern struct { - double chi[ReaxParams::nsort]; - double eta[ReaxParams::nsort]; - double gam[ReaxParams::nsort]; -} FORTRAN(cbkchb,CBKCHB); diff --git a/src/REAX/reax_cbkd.h b/src/REAX/reax_cbkd.h deleted file mode 100644 index 78536665e4..0000000000 --- a/src/REAX/reax_cbkd.h +++ /dev/null @@ -1,12 +0,0 @@ -/////:EOH~ -extern struct { - double d[3*ReaxParams::nat]; double estrain[ReaxParams::nat]; -} FORTRAN(cbkd,CBKD); - -extern struct { - double atomvirial[6*ReaxParams::nat]; - double virial[6]; - int Lvirial; - int Latomvirial; -} FORTRAN(cbkvirial,CBKVIRIAL); - diff --git a/src/REAX/reax_cbkia.h b/src/REAX/reax_cbkia.h deleted file mode 100644 index 3681005292..0000000000 --- a/src/REAX/reax_cbkia.h +++ /dev/null @@ -1,6 +0,0 @@ -/////:EOH~ -extern struct { - int ia[ReaxParams::nat*(ReaxParams::mbond+3)]; - int iag[ReaxParams::nat*(ReaxParams::mbond+3)]; -} FORTRAN(cbkia,CBKIA); - diff --git a/src/REAX/reax_cbklonpar.h b/src/REAX/reax_cbklonpar.h deleted file mode 100644 index 4bf27522e9..0000000000 --- a/src/REAX/reax_cbklonpar.h +++ /dev/null @@ -1,5 +0,0 @@ -/////:EOH~ -extern struct { - double vlp[ReaxParams::nat]; - double dvlpdsbo[ReaxParams::nat]; -} FORTRAN(cbklonpar,CBKLONPAR); diff --git a/src/REAX/reax_cbknubon2.h b/src/REAX/reax_cbknubon2.h deleted file mode 100644 index 22ab3d569a..0000000000 --- a/src/REAX/reax_cbknubon2.h +++ /dev/null @@ -1,6 +0,0 @@ -/////:EOH~ -extern struct { - int nubon1[ReaxParams::nat*(ReaxParams::mbond)]; - int nubon2[ReaxParams::nat*(ReaxParams::mbond)]; -} FORTRAN(cbknubon2,CBKNUBON2); - diff --git a/src/REAX/reax_cbkpairs.h b/src/REAX/reax_cbkpairs.h deleted file mode 100644 index 3f7653dd89..0000000000 --- a/src/REAX/reax_cbkpairs.h +++ /dev/null @@ -1,15 +0,0 @@ -/////:EOH~ -extern struct { - int nvl1[ReaxParams::nneighmax * ReaxParams::nat]; - int nvl2[ReaxParams::nneighmax * ReaxParams::nat]; - int nvpair; - int nvlself; -} FORTRAN(cbkpairs,CBKPAIRS); - -extern struct { - int nvlbo[ReaxParams::nneighmax * ReaxParams::nat]; -} FORTRAN(cbknvlbo,CBKNVLBO); - -extern struct { - int nvlown[ReaxParams::nneighmax * ReaxParams::nat]; -} FORTRAN(cbknvlown,CBKNVLOWN); diff --git a/src/REAX/reax_cbkqa.h b/src/REAX/reax_cbkqa.h deleted file mode 100644 index cc54a8558f..0000000000 --- a/src/REAX/reax_cbkqa.h +++ /dev/null @@ -1,5 +0,0 @@ -/////:EOH~ -extern struct { - char qa[20*ReaxParams::nattot+10]; -} FORTRAN(cbkqa,CBKQA); - diff --git a/src/REAX/reax_energies.h b/src/REAX/reax_energies.h deleted file mode 100644 index 76ad6990da..0000000000 --- a/src/REAX/reax_energies.h +++ /dev/null @@ -1,25 +0,0 @@ -/////:EOH~ - -extern struct { - double eb; - double eoop; - double epen; - double estrc; - double deda[3]; - double pressu; - double efi; - double elp; - double emol; - double ea; - double eres; - double et; - double eradbo; - double ev; - double eco; - double ecoa; - double ehb; - double sw; - double ew; - double ep; - double ekin; -} FORTRAN(cbkenergies,CBKENERGIES); diff --git a/src/REAX/reax_fortran.h b/src/REAX/reax_fortran.h deleted file mode 100644 index c79ae5ab93..0000000000 --- a/src/REAX/reax_fortran.h +++ /dev/null @@ -1,31 +0,0 @@ -/////:EOH~ - -/* CONS(a,b) should return ab, the concatenation - of its arguments */ -#if __STDC__ -#define CONS(a,b) a##b -#else -#define CONS(a,b) a/**/b -#endif - -#ifdef _ARDENT -#define FORTRAN(lcname,ucname) ucname -#endif - -#ifdef _IBM -#define FORTRAN(lcname,ucname) lcname -#endif - -#ifdef _HP -#define FORTRAN(lcname,ucname) lcname -#endif - -#ifdef _F2C_LINUX -#define FORTRAN(lcname,ucname) CONS(lcname,__) -#endif - -#ifndef FORTRAN -#define FORTRAN(lcname,ucname) CONS(lcname,_) -#endif - - diff --git a/src/REAX/reax_functions.h b/src/REAX/reax_functions.h deleted file mode 100644 index 5a0a3acfb5..0000000000 --- a/src/REAX/reax_functions.h +++ /dev/null @@ -1,37 +0,0 @@ -/////:EOH~ - -extern "C" void FORTRAN(readc,READC)(); -extern "C" void FORTRAN(init,INIT)(); -extern "C" void FORTRAN(ffinpt,FFINPT)(); -extern "C" void FORTRAN(tap7th,TAP7TH)(); -extern "C" void FORTRAN(taper,TAPER)(double*,double*); -extern "C" void FORTRAN(readgeo,READGEO)(); -extern "C" void FORTRAN(srtatom,SRTATOM)(); -extern "C" void FORTRAN(vlist,VLIST) (); -extern "C" void FORTRAN(srtbon1,SRTBON1)(int*,int*,double*); -extern "C" void FORTRAN(molec,MOLEC)(); -extern "C" void FORTRAN(encalc,ENCALC)(); -extern "C" void FORTRAN(getswb,GETSWB)(double*); -extern "C" void FORTRAN(getswa,GETSWA)(double*); -extern "C" void FORTRAN(getvrange,GET_VRANGE)(double*); -extern "C" void FORTRAN(getnvlist,GET_NVLIST)(int*); -extern "C" void FORTRAN(getvlbora,GETVLBORA)(double*); -extern "C" void FORTRAN(cgsolve,CGSOLVE)(int*,double*,int*,double*,double*,int*); -extern "C" void FORTRAN(getnval,GETNVAL)(int*); -extern "C" void FORTRAN(getntor,GETNTOR)(int*); -extern "C" void FORTRAN(getnhb,GETNHB)(int*); -extern "C" void FORTRAN(getnbonall,GETNBONALL)(int*); -extern "C" void FORTRAN(getnneighmax,GETNNEIGHMAX)(int*); -extern "C" void FORTRAN(getnat,GETNAT)(int*); -extern "C" void FORTRAN(getnattot,GETNATTOT)(int*); -extern "C" void FORTRAN(getnsort,GETNSORT)(int*); -extern "C" void FORTRAN(getmbond,GETMBOND)(int*); -extern "C" void FORTRAN(getnso,GETNSO)(int*); -extern "C" void FORTRAN(setngeofor,SETNGEOFOR)(int*); -extern "C" void FORTRAN(mdsav,MDSAV)(int*); -extern "C" void FORTRAN(getnsbmax,GETNSBMAX)(int*); -extern "C" void FORTRAN(getnsbma2,GETNSBMA2)(int*); -extern "C" void FORTRAN(getcutof3,GETCUTOF3)(double*); - - - diff --git a/src/REAX/reax_params.h b/src/REAX/reax_params.h deleted file mode 100644 index 7e57a24ce4..0000000000 --- a/src/REAX/reax_params.h +++ /dev/null @@ -1,14 +0,0 @@ -/////:EOH~ - -#include "reax_defs.h" - -class ReaxParams { - public: - enum {nneighmax=NNEIGHMAXDEF, - nat=NATDEF, - nattot=NATTOTDEF, - nsort=NSORTDEF, - mbond=MBONDDEF, - nbomax=NBOMAXDEF, - }; -}; diff --git a/src/REAX/reax_small.h b/src/REAX/reax_small.h deleted file mode 100644 index 1cec2473ec..0000000000 --- a/src/REAX/reax_small.h +++ /dev/null @@ -1,15 +0,0 @@ -/////:EOH~ - -extern struct { - double tset; - double dseed; - double tempmd; - double ts2; - double ts22; - int nmolo; - int nmolo5; - int nbon; - int na; - int namov; - int na_local;} FORTRAN(rsmall,RSMALL); -