diff -Naur lammps-6Jun09/doc/pair_reax.txt lammps-22Jun09/doc/pair_reax.txt --- lammps-6Jun09/doc/pair_reax.txt 2009-06-02 11:52:25.000000000 -0600 +++ lammps-22Jun09/doc/pair_reax.txt 2009-06-09 13:21:42.000000000 -0600 @@ -36,7 +36,8 @@ LAMMPS provides a ReaxFF potential file in its potentials dir, namely potentials/ffield.reax. Its format is identical to that used by van Duin and co-workers. It contains parameterizations for the following -elements: C, H, O, N, S, and Si. You can use your another file in +elements: C, H, O, N, S (Si has been temporarily removed). +You can use another file in place of it, and ReaxFF files with parameterizations for other elements or for specific chemical systems may be available elsewhere. diff -Naur lammps-6Jun09/examples/reax/ffield.reax lammps-22Jun09/examples/reax/ffield.reax --- lammps-6Jun09/examples/reax/ffield.reax 2009-06-01 14:16:04.000000000 -0600 +++ lammps-22Jun09/examples/reax/ffield.reax 2009-06-09 13:19:29.000000000 -0600 @@ -39,7 +39,7 @@ 5.0000 !Molecular energy (not used) 0.0000 !Molecular energy (not used) 3.6942 !Valency angle conjugation parameter - 6 ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;# + 5 ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;# alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u. cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u. ov/un;val1;n.u.;val3,vval4 @@ -63,11 +63,7 @@ 9.9575 4.9055 4.0000 52.9998 112.1416 6.5000 8.2545 2.0000 1.4601 9.7177 -2.3700 5.7487 23.2859 12.7147 0.9745 0.0000 -11.0000 2.7466 1.0338 4.0000 2.8793 0.0000 0.0000 0.0000 - Si 2.0276 4.0000 28.0600 2.2042 0.1322 0.8218 1.5758 4.0000 - 11.9413 2.0618 4.0000 11.8211 136.4845 1.8038 7.3852 0.0000 - -1.0000 0.0000 -2.3700 6.4918 8.5961 0.2368 0.8563 0.0000 - -3.8112 3.1873 1.0338 4.0000 2.5791 0.0000 0.0000 0.0000 - 18 ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6 + 15 ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6 pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr 1 1 145.4070 103.0681 73.7841 0.2176 -0.7816 1.0000 28.4167 0.3217 0.1111 -0.1940 8.6733 1.0000 -0.0994 5.9724 1.0000 0.0000 @@ -99,22 +95,14 @@ 0.3296 -0.3153 9.1227 1.0000 -0.1805 5.6864 1.0000 0.0000 5 5 96.1871 93.7006 68.6860 0.0955 -0.4781 1.0000 17.8574 0.6000 0.2723 -0.2373 9.7875 1.0000 -0.0950 6.4757 1.0000 0.0000 - 6 6 109.1904 70.8314 30.0000 0.2765 -0.3000 1.0000 16.0000 0.1583 - 0.2804 -0.1994 8.1117 1.0000 -0.0675 8.2993 0.0000 0.0000 - 2 6 137.1002 0.0000 0.0000 -0.1902 0.0000 1.0000 6.0000 0.4256 - 17.7186 1.0000 0.0000 1.0000 -0.0377 6.4281 0.0000 0.0000 - 3 6 191.1743 52.0733 43.3991 -0.2584 -0.3000 1.0000 36.0000 0.8764 - 1.0248 -0.3658 4.2151 1.0000 -0.5004 4.2605 1.0000 0.0000 - 8 ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2 + 6 ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2 1 2 0.0455 1.7218 10.4236 1.0379 -1.0000 -1.0000 2 3 0.0469 1.9185 10.3707 0.9406 -1.0000 -1.0000 2 4 0.0999 1.8372 9.6539 0.9692 -1.0000 -1.0000 1 3 0.1186 1.9820 9.5927 1.2936 1.1203 1.0805 1 4 0.1486 1.8922 9.7989 1.3746 1.2091 1.1427 3 4 0.1051 2.0060 10.0691 1.3307 1.1034 1.0060 - 2 6 0.0470 1.6738 11.6877 1.1931 -1.0000 -1.0000 - 3 6 0.1263 1.8163 10.6833 1.6266 1.2052 -1.0000 - 62 ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2 + 50 ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2 1 1 1 70.0265 13.6338 2.1884 0.0000 0.1676 26.3587 1.0400 1 1 2 69.7786 10.3544 8.4326 0.0000 0.1153 0.0000 1.0400 2 1 2 74.6020 11.8629 2.9294 0.0000 0.1367 0.0000 1.0400 @@ -164,20 +152,8 @@ 1 5 5 85.3644 36.9951 2.0903 0.1463 0.0559 0.0000 1.0400 2 5 2 93.1959 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400 2 5 5 84.3331 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400 - 6 6 6 69.3456 21.7361 1.4283 0.0000 -0.2101 0.0000 1.3241 - 2 6 6 75.6168 21.5317 1.0435 0.0000 2.5179 0.0000 1.0400 - 2 6 2 78.3939 20.9772 0.8630 0.0000 2.8421 0.0000 1.0400 - 3 6 6 70.3016 15.4081 1.3267 0.0000 2.1459 0.0000 1.0400 - 2 6 3 73.8232 16.6592 3.7425 0.0000 0.8613 0.0000 1.0400 - 3 6 3 90.0344 7.7656 1.7264 0.0000 0.7689 0.0000 1.0400 - 6 3 6 22.1715 3.6615 0.3160 0.0000 4.1125 0.0000 1.0400 - 2 3 6 83.7634 5.6693 2.7780 0.0000 1.6982 0.0000 1.0400 - 3 3 6 73.4663 25.0761 0.9143 0.0000 2.2466 0.0000 1.0400 - 2 2 6 0.0000 47.1300 6.0000 0.0000 1.6371 0.0000 1.0400 - 6 2 6 0.0000 31.5209 6.0000 0.0000 1.6371 0.0000 1.0400 - 3 2 6 0.0000 31.0427 4.5625 0.0000 1.6371 0.0000 1.0400 2 2 5 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400 - 20 ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n + 17 ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n 1 1 1 1 0.0000 23.2168 0.1811 -4.6220 -1.9387 0.0000 0.0000 1 1 1 2 0.0000 45.7984 0.3590 -5.7106 -2.9459 0.0000 0.0000 2 1 1 2 0.0000 44.6445 0.3486 -5.1725 -0.8717 0.0000 0.0000 @@ -195,9 +171,6 @@ 0 1 5 0 3.3423 30.3435 0.0365 -2.7171 0.0000 0.0000 0.0000 0 5 5 0 -0.0555 -42.7738 0.1515 -2.2056 0.0000 0.0000 0.0000 0 2 5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 - 0 6 6 0 0.0000 0.0000 0.1200 -2.4426 0.0000 0.0000 0.0000 - 0 2 6 0 0.0000 0.0000 0.1200 -2.4847 0.0000 0.0000 0.0000 - 0 3 6 0 0.0000 0.0000 0.1200 -2.4703 0.0000 0.0000 0.0000 9 ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1 3 2 3 2.0431 -6.6813 3.5000 1.7295 3 2 4 1.6740 -10.9581 3.5000 1.7295 diff -Naur lammps-6Jun09/src/MOLECULE/atom_vec_angle.cpp lammps-22Jun09/src/MOLECULE/atom_vec_angle.cpp --- lammps-6Jun09/src/MOLECULE/atom_vec_angle.cpp 2007-10-04 11:57:04.000000000 -0600 +++ lammps-22Jun09/src/MOLECULE/atom_vec_angle.cpp 2009-06-22 13:59:24.000000000 -0600 @@ -288,7 +288,7 @@ int AtomVecAngle::pack_border_one(int i, double *buf) { - buf[0] = v[i][0]; + buf[0] = molecule[i]; return 1; } diff -Naur lammps-6Jun09/src/REAX/pair_reax.cpp lammps-22Jun09/src/REAX/pair_reax.cpp --- lammps-6Jun09/src/REAX/pair_reax.cpp 2009-06-02 12:09:07.000000000 -0600 +++ lammps-22Jun09/src/REAX/pair_reax.cpp 2009-06-12 11:57:10.000000000 -0600 @@ -1030,3 +1030,62 @@ bytes += matmax * sizeof(double); return bytes; } + +/* ---------------------------------------------------------------------- + print out the itemized energies to the log file + this is not currently called, but should + be made available to the user in some way +------------------------------------------------------------------------- */ + +void PairREAX::output_itemized_energy(double energy_charge_equilibration) +{ + double etmp[14],etmp2[14]; + + etmp[0] = FORTRAN(cbkenergies, CBKENERGIES).eb; + etmp[1] = FORTRAN(cbkenergies, CBKENERGIES).ea; + etmp[2] = FORTRAN(cbkenergies, CBKENERGIES).elp; + etmp[3] = FORTRAN(cbkenergies, CBKENERGIES).emol; + etmp[4] = FORTRAN(cbkenergies, CBKENERGIES).ev; + etmp[5] = FORTRAN(cbkenergies, CBKENERGIES).epen; + etmp[6] = FORTRAN(cbkenergies, CBKENERGIES).ecoa; + etmp[7] = FORTRAN(cbkenergies, CBKENERGIES).ehb; + etmp[8] = FORTRAN(cbkenergies, CBKENERGIES).et; + etmp[9] = FORTRAN(cbkenergies, CBKENERGIES).eco; + etmp[10] = FORTRAN(cbkenergies, CBKENERGIES).ew; + etmp[11] = FORTRAN(cbkenergies, CBKENERGIES).ep; + etmp[12] = FORTRAN(cbkenergies, CBKENERGIES).efi; + etmp[13] = energy_charge_equilibration; + + MPI_Allreduce(etmp,etmp2,14,MPI_DOUBLE,MPI_SUM,world); + + if (comm->me == 0) { + fprintf(screen,"eb = %g \n",etmp2[0]); + fprintf(screen,"ea = %g \n",etmp2[1]); + fprintf(screen,"elp = %g \n",etmp2[2]); + fprintf(screen,"emol = %g \n",etmp2[3]); + fprintf(screen,"ecoa = %g \n",etmp2[4]); + fprintf(screen,"epen = %g \n",etmp2[5]); + fprintf(screen,"ecoa = %g \n",etmp2[6]); + fprintf(screen,"ehb = %g \n",etmp2[7]); + fprintf(screen,"et = %g \n",etmp2[8]); + fprintf(screen,"eco = %g \n",etmp2[9]); + fprintf(screen,"ew = %g \n",etmp2[10]); + fprintf(screen,"ep = %g \n",etmp2[11]); + fprintf(screen,"efi = %g \n",etmp2[12]); + fprintf(screen,"eqeq = %g \n",etmp2[13]); + fprintf(logfile,"eb = %g \n",etmp2[0]); + fprintf(logfile,"ea = %g \n",etmp2[1]); + fprintf(logfile,"elp = %g \n",etmp2[2]); + fprintf(logfile,"emol = %g \n",etmp2[3]); + fprintf(logfile,"ecoa = %g \n",etmp2[4]); + fprintf(logfile,"epen = %g \n",etmp2[5]); + fprintf(logfile,"ecoa = %g \n",etmp2[6]); + fprintf(logfile,"ehb = %g \n",etmp2[7]); + fprintf(logfile,"et = %g \n",etmp2[8]); + fprintf(logfile,"eco = %g \n",etmp2[9]); + fprintf(logfile,"ew = %g \n",etmp2[10]); + fprintf(logfile,"ep = %g \n",etmp2[11]); + fprintf(logfile,"efi = %g \n",etmp2[12]); + fprintf(logfile,"eqeq = %g \n",etmp2[13]); + } +} diff -Naur lammps-6Jun09/src/REAX/pair_reax.h lammps-22Jun09/src/REAX/pair_reax.h --- lammps-6Jun09/src/REAX/pair_reax.h 2009-02-12 13:45:31.000000000 -0700 +++ lammps-22Jun09/src/REAX/pair_reax.h 2009-06-12 11:57:10.000000000 -0600 @@ -82,6 +82,7 @@ int[], double[], double[]); void charge_reax(const int &, const int &, double[], double[], int[], int[], double[]); + void output_itemized_energy(double); }; } diff -Naur lammps-6Jun09/src/pair.cpp lammps-22Jun09/src/pair.cpp --- lammps-6Jun09/src/pair.cpp 2009-02-12 13:42:48.000000000 -0700 +++ lammps-22Jun09/src/pair.cpp 2009-06-08 12:29:00.000000000 -0600 @@ -691,6 +691,72 @@ } /* ---------------------------------------------------------------------- + tally virial into global and per-atom accumulators + called by pair lubricate potential with 6 tensor components +------------------------------------------------------------------------- */ + +void Pair::v_tally_tensor(int i, int j, int nlocal, int newton_pair, + double vxx, double vyy, double vzz, + double vxy, double vxz, double vyz) +{ + double v[6]; + + v[0] = vxx; + v[1] = vyy; + v[2] = vzz; + v[3] = vxy; + v[4] = vxz; + v[5] = vyz; + + if (vflag_global) { + if (newton_pair) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + if (j < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_pair || i < nlocal) { + vatom[i][0] += 0.5*v[0]; + vatom[i][1] += 0.5*v[1]; + vatom[i][2] += 0.5*v[2]; + vatom[i][3] += 0.5*v[3]; + vatom[i][4] += 0.5*v[4]; + vatom[i][5] += 0.5*v[5]; + } + if (newton_pair || j < nlocal) { + vatom[j][0] += 0.5*v[0]; + vatom[j][1] += 0.5*v[1]; + vatom[j][2] += 0.5*v[2]; + vatom[j][3] += 0.5*v[3]; + vatom[j][4] += 0.5*v[4]; + vatom[j][5] += 0.5*v[5]; + } + } +} + +/* ---------------------------------------------------------------------- compute global pair virial via summing F dot r over own & ghost atoms at this point, only pairwise forces have been accumulated in atom->f ------------------------------------------------------------------------- */ diff -Naur lammps-6Jun09/src/pair.h lammps-22Jun09/src/pair.h --- lammps-6Jun09/src/pair.h 2009-01-05 15:26:08.000000000 -0700 +++ lammps-22Jun09/src/pair.h 2009-06-08 12:29:00.000000000 -0600 @@ -127,6 +127,8 @@ void v_tally3(int, int, int, double *, double *, double *, double *); void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *); + void v_tally_tensor(int, int, int, int, + double, double, double, double, double, double); void virial_compute(); };