diff --git a/doc/src/pair_dipole.txt b/doc/src/pair_dipole.txt index 14554f1df..1c557fbc6 100644 --- a/doc/src/pair_dipole.txt +++ b/doc/src/pair_dipole.txt @@ -1,263 +1,265 @@ "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c :link(lws,http://lammps.sandia.gov) :link(ld,Manual.html) :link(lc,Section_commands.html#comm) :line pair_style lj/cut/dipole/cut command :h3 pair_style lj/cut/dipole/cut/gpu command :h3 pair_style lj/cut/dipole/cut/omp command :h3 pair_style lj/sf/dipole/sf command :h3 pair_style lj/sf/dipole/sf/gpu command :h3 pair_style lj/sf/dipole/sf/omp command :h3 pair_style lj/cut/dipole/long command :h3 pair_style lj/long/dipole/long command :h3 [Syntax:] pair_style lj/cut/dipole/cut cutoff (cutoff2) pair_style lj/sf/dipole/sf cutoff (cutoff2) pair_style lj/cut/dipole/long cutoff (cutoff2) pair_style lj/long/dipole/long flag_lj flag_coul cutoff (cutoff2) :pre cutoff = global cutoff LJ (and Coulombic if only 1 arg) (distance units) :ulb,l cutoff2 = global cutoff for Coulombic and dipole (optional) (distance units) :l flag_lj = {long} or {cut} or {off} :l {long} = use long-range damping on dispersion 1/r^6 term {cut} = use a cutoff on dispersion 1/r^6 term {off} = omit disperion 1/r^6 term entirely :pre flag_coul = {long} or {off} :l {long} = use long-range damping on Coulombic 1/r and point-dipole terms {off} = omit Coulombic and point-dipole terms entirely :pre :ule [Examples:] pair_style lj/cut/dipole/cut 10.0 pair_coeff * * 1.0 1.0 pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre -pair_style lj/sf/dipole/sf 9.0 +pair_style lj/sf/dipole/sf 3.0 3.0 pair_coeff * * 1.0 1.0 -pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre +pair_coeff * * 1.0 1.0 scale 0.8 +pair_coeff 2 3 1.0 1.0 2.5 4.0 +pair_coeff * * 1.0 1.0 2.5 4.0 scale 0.8 :pre pair_style lj/cut/dipole/long 10.0 pair_coeff * * 1.0 1.0 pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre pair_style lj/long/dipole/long long long 3.5 10.0 pair_coeff * * 1.0 1.0 pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre [Description:] Style {lj/cut/dipole/cut} computes interactions between pairs of particles that each have a charge and/or a point dipole moment. In addition to the usual Lennard-Jones interaction between the particles (Elj) the charge-charge (Eqq), charge-dipole (Eqp), and dipole-dipole (Epp) interactions are computed by these formulas for the energy (E), force (F), and torque (T) between particles I and J. :c,image(Eqs/pair_dipole.jpg) where qi and qj are the charges on the two particles, pi and pj are the dipole moment vectors of the two particles, r is their separation distance, and the vector r = Ri - Rj is the separation vector between the two particles. Note that Eqq and Fqq are simply Coulombic energy and force, Fij = -Fji as symmetric forces, and Tij != -Tji since the torques do not act symmetrically. These formulas are discussed in "(Allen)"_#Allen and in "(Toukmaji)"_#Toukmaji. Style {lj/sf/dipole/sf} computes "shifted-force" interactions between pairs of particles that each have a charge and/or a point dipole moment. In general, a shifted-force potential is a (sligthly) modified potential containing extra terms that make both the energy and its derivative go to zero at the cutoff distance; this removes (cutoff-related) problems in energy conservation and any numerical instability in the equations of motion "(Allen)"_#Allen. Shifted-force interactions for the Lennard-Jones (E_LJ), charge-charge (Eqq), charge-dipole (Eqp), dipole-charge (Epq) and dipole-dipole (Epp) potentials are computed by these formulas for the energy (E), force (F), and torque (T) between particles I and J: :c,image(Eqs/pair_dipole_sf.jpg) :c,image(Eqs/pair_dipole_sf2.jpg) where epsilon and sigma are the standard LJ parameters, r_c is the cutoff, qi and qj are the charges on the two particles, pi and pj are the dipole moment vectors of the two particles, r is their separation distance, and the vector r = Ri - Rj is the separation vector between the two particles. Note that Eqq and Fqq are simply Coulombic energy and force, Fij = -Fji as symmetric forces, and Tij != -Tji since the torques do not act symmetrically. The shifted-force formula for the Lennard-Jones potential is reported in "(Stoddard)"_#Stoddard. The original (unshifted) formulas for the electrostatic potentials, forces and torques can be found in "(Price)"_#Price. The shifted-force electrostatic potentials have been obtained by applying equation 5.13 of "(Allen)"_#Allen. The formulas for the corresponding forces and torques have been obtained by applying the 'chain rule' as in appendix C.3 of "(Allen)"_#Allen. If one cutoff is specified in the pair_style command, it is used for both the LJ and Coulombic (q,p) terms. If two cutoffs are specified, they are used as cutoffs for the LJ and Coulombic (q,p) terms respectively. Style {lj/cut/dipole/long} computes long-range point-dipole interactions as discussed in "(Toukmaji)"_#Toukmaji. Dipole-dipole, dipole-charge, and charge-charge interactions are all supported, along with the standard 12/6 Lennard-Jones interactions, which are computed with a cutoff. A "kspace_style"_kspace_style.html must be defined to use this pair style. Currently, only "kspace_style ewald/disp"_kspace_style.html support long-range point-dipole interactions. Style {lj/long/dipole/long} also computes point-dipole interactions as discussed in "(Toukmaji)"_#Toukmaji. Long-range dipole-dipole, dipole-charge, and charge-charge interactions are all supported, along with the standard 12/6 Lennard-Jones interactions. LJ interactions can be cutoff or long-ranged. For style {lj/long/dipole/long}, if {flag_lj} is set to {long}, no cutoff is used on the LJ 1/r^6 dispersion term. The long-range portion is calculated by using the "kspace_style ewald_disp"_kspace_style.html command. The specified LJ cutoff then determines which portion of the LJ interactions are computed directly by the pair potential versus which part is computed in reciprocal space via the Kspace style. If {flag_lj} is set to {cut}, the LJ interactions are simply cutoff, as with "pair_style lj/cut"_pair_lj.html. If {flag_lj} is set to {off}, LJ interactions are not computed at all. If {flag_coul} is set to {long}, no cutoff is used on the Coulombic or dipole interactions. The long-range portion is calculated by using {ewald_disp} of the "kspace_style"_kspace_style.html command. If {flag_coul} is set to {off}, Coulombic and dipole interactions are not computed at all. Atoms with dipole moments should be integrated using the "fix nve/sphere update dipole"_fix_nve_sphere.html command to rotate the dipole moments. The {omega} option on the "fix langevin"_fix_langevin.html command can be used to thermostat the rotational motion. The "compute temp/sphere"_compute_temp_sphere.html command can be used to monitor the temperature, since it includes rotational degrees of freedom. The "atom_style dipole"_atom_style.html command should be used since it defines the point dipoles and their rotational state. The magnitude of the dipole moment for each type of particle can be defined by the "dipole"_dipole.html command or in the "Dipoles" section of the data file read in by the "read_data"_read_data.html command. Their initial orientation can be defined by the "set dipole"_set.html command or in the "Atoms" section of the data file. The following coefficients must be defined for each pair of atoms types via the "pair_coeff"_pair_coeff.html command as in the examples above, or in the data file or restart files read by the "read_data"_read_data.html or "read_restart"_read_restart.html commands, or by mixing as described below: epsilon (energy units) sigma (distance units) cutoff1 (distance units) cutoff2 (distance units) :ul The latter 2 coefficients are optional. If not specified, the global LJ and Coulombic cutoffs specified in the pair_style command are used. If only one cutoff is specified, it is used as the cutoff for both LJ and Coulombic interactions for this type pair. If both coefficients are specified, they are used as the LJ and Coulombic cutoffs for this type pair. :line Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are functionally the same as the corresponding style without the suffix. They have been optimized to run faster, depending on your available hardware, as discussed in "Section_accelerate"_Section_accelerate.html of the manual. The accelerated styles take the same arguments and should produce the same results, except for round-off and precision issues. These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, USER-OMP and OPT packages, respectively. They are only enabled if LAMMPS was built with those packages. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. You can specify the accelerated styles explicitly in your input script by including their suffix, or you can use the "-suffix command-line switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can use the "suffix"_suffix.html command in your input script. See "Section_accelerate"_Section_accelerate.html of the manual for more instructions on how to use the accelerated styles effectively. :line [Mixing, shift, table, tail correction, restart, rRESPA info]: For atom type pairs I,J and I != J, the epsilon and sigma coefficients and cutoff distances for this pair style can be mixed. The default mix value is {geometric}. See the "pair_modify" command for details. For atom type pairs I,J and I != J, the A, sigma, d1, and d2 coefficients and cutoff distance for this pair style can be mixed. A is an energy value mixed like a LJ epsilon. D1 and d2 are distance values and are mixed like sigma. The default mix value is {geometric}. See the "pair_modify" command for details. This pair style does not support the "pair_modify"_pair_modify.html shift option for the energy of the Lennard-Jones portion of the pair interaction; such energy goes to zero at the cutoff by construction. The "pair_modify"_pair_modify.html table option is not relevant for this pair style. This pair style does not support the "pair_modify"_pair_modify.html tail option for adding long-range tail corrections to energy and pressure. This pair style writes its information to "binary restart files"_restart.html, so pair_style and pair_coeff commands do not need to be specified in an input script that reads a restart file. This pair style can only be used via the {pair} keyword of the "run_style respa"_run_style.html command. It does not support the {inner}, {middle}, {outer} keywords. [Restrictions:] The {lj/cut/dipole/cut}, {lj/cut/dipole/long}, and {lj/long/dipole/long} styles are part of the DIPOLE package. They are only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. The {lj/sf/dipole/sf} style is part of the USER-MISC package. It is only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. Using dipole pair styles with {electron} "units"_units.html is not currently supported. [Related commands:] "pair_coeff"_pair_coeff.html [Default:] none :line :link(Allen) [(Allen)] Allen and Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. :link(Toukmaji) [(Toukmaji)] Toukmaji, Sagui, Board, and Darden, J Chem Phys, 113, 10913 (2000). :link(Stoddard) [(Stoddard)] Stoddard and Ford, Phys Rev A, 8, 1504 (1973). :link(Price) [(Price)] Price, Stone and Alderton, Mol Phys, 52, 987 (1984). diff --git a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp index 072a6c994..b88f72bac 100644 --- a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp +++ b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp @@ -1,621 +1,642 @@ /* ---------------------------------------------------------------------- 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 authors: Mario Orsi (QMUL), m.orsi@qmul.ac.uk Samuel Genheden (University of Southampton) ------------------------------------------------------------------------- */ #include #include #include "pair_lj_sf_dipole_sf.h" #include "atom.h" #include "neighbor.h" #include "neigh_list.h" #include "comm.h" #include "force.h" #include "memory.h" #include "error.h" #include "update.h" #include using namespace LAMMPS_NS; +static int warn_single = 0; + /* ---------------------------------------------------------------------- */ PairLJSFDipoleSF::PairLJSFDipoleSF(LAMMPS *lmp) : Pair(lmp) { } /* ---------------------------------------------------------------------- */ PairLJSFDipoleSF::~PairLJSFDipoleSF() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut_lj); memory->destroy(cut_ljsq); memory->destroy(cut_coul); memory->destroy(cut_coulsq); memory->destroy(epsilon); memory->destroy(sigma); memory->destroy(lj1); memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); memory->destroy(scale); } } /* ---------------------------------------------------------------------- */ void PairLJSFDipoleSF::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fx,fy,fz; double rsq,rinv,r2inv,r6inv,r3inv,r5inv; double forcecoulx,forcecouly,forcecoulz,crossx,crossy,crossz; double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul; double fq,pdotp,pidotr,pjdotr,pre1,pre2,pre3,pre4; double forcelj,factor_coul,factor_lj; double presf,afac,bfac,pqfac,qpfac,forceljcut,forceljsf; double aforcecoulx,aforcecouly,aforcecoulz; double bforcecoulx,bforcecouly,bforcecoulz; double rcutlj2inv, rcutcoul2inv,rcutlj6inv; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; double *q = atom->q; double **mu = atom->mu; double **torque = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; -// The global scaling parameters aren't used anymore double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; factor_coul = special_coul[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; rinv = sqrt(r2inv); // atom can have both a charge and dipole // i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole forcecoulx = forcecouly = forcecoulz = 0.0; tixcoul = tiycoul = tizcoul = 0.0; tjxcoul = tjycoul = tjzcoul = 0.0; if (rsq < cut_coulsq[itype][jtype]) { rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; if (qtmp != 0.0 && q[j] != 0.0) { pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv); forcecoulx += pre1*delx; forcecouly += pre1*dely; forcecoulz += pre1*delz; } if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; afac = 1.0 - rsq*rsq * rcutcoul2inv*rcutcoul2inv; pre1 = afac * ( pdotp - 3.0 * r2inv * pidotr * pjdotr ); aforcecoulx = pre1*delx; aforcecouly = pre1*dely; aforcecoulz = pre1*delz; bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv + 3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv; presf = 2.0 * r2inv * pidotr * pjdotr; bforcecoulx = bfac * (pjdotr*mu[i][0]+pidotr*mu[j][0]-presf*delx); bforcecouly = bfac * (pjdotr*mu[i][1]+pidotr*mu[j][1]-presf*dely); bforcecoulz = bfac * (pjdotr*mu[i][2]+pidotr*mu[j][2]-presf*delz); forcecoulx += 3.0 * r5inv * ( aforcecoulx + bforcecoulx ); forcecouly += 3.0 * r5inv * ( aforcecouly + bforcecouly ); forcecoulz += 3.0 * r5inv * ( aforcecoulz + bforcecoulz ); pre2 = 3.0 * bfac * r5inv * pjdotr; pre3 = 3.0 * bfac * r5inv * pidotr; pre4 = -bfac * r3inv; crossx = pre4 * (mu[i][1]*mu[j][2] - mu[i][2]*mu[j][1]); crossy = pre4 * (mu[i][2]*mu[j][0] - mu[i][0]*mu[j][2]); crossz = pre4 * (mu[i][0]*mu[j][1] - mu[i][1]*mu[j][0]); tixcoul += crossx + pre2 * (mu[i][1]*delz - mu[i][2]*dely); tiycoul += crossy + pre2 * (mu[i][2]*delx - mu[i][0]*delz); tizcoul += crossz + pre2 * (mu[i][0]*dely - mu[i][1]*delx); tjxcoul += -crossx + pre3 * (mu[j][1]*delz - mu[j][2]*dely); tjycoul += -crossy + pre3 * (mu[j][2]*delx - mu[j][0]*delz); tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx); } if (mu[i][3] > 0.0 && q[j] != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv); pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = q[j] * r3inv * pqfac; forcecoulx += pre2*mu[i][0] - pre1*delx; forcecouly += pre2*mu[i][1] - pre1*dely; forcecoulz += pre2*mu[i][2] - pre1*delz; tixcoul += pre2 * (mu[i][1]*delz - mu[i][2]*dely); tiycoul += pre2 * (mu[i][2]*delx - mu[i][0]*delz); tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx); } if (mu[j][3] > 0.0 && qtmp != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv); qpfac = 1.0 - 3.0*rsq*rcutcoul2inv + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = qtmp * r3inv * qpfac; forcecoulx += pre1*delx - pre2*mu[j][0]; forcecouly += pre1*dely - pre2*mu[j][1]; forcecoulz += pre1*delz - pre2*mu[j][2]; tjxcoul += -pre2 * (mu[j][1]*delz - mu[j][2]*dely); tjycoul += -pre2 * (mu[j][2]*delx - mu[j][0]*delz); tjzcoul += -pre2 * (mu[j][0]*dely - mu[j][1]*delx); } } // LJ interaction if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; forceljcut = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype])*r2inv; rcutlj2inv = 1.0 / cut_ljsq[itype][jtype]; rcutlj6inv = rcutlj2inv * rcutlj2inv * rcutlj2inv; forceljsf = (lj1[itype][jtype]*rcutlj6inv - lj2[itype][jtype]) * rcutlj6inv * rcutlj2inv; forcelj = factor_lj * (forceljcut - forceljsf); } else forcelj = 0.0; // total force fq = factor_coul*qqrd2e*scale[itype][jtype]; fx = fq*forcecoulx + delx*forcelj; fy = fq*forcecouly + dely*forcelj; fz = fq*forcecoulz + delz*forcelj; // force & torque accumulation f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; torque[i][0] += fq*tixcoul; torque[i][1] += fq*tiycoul; torque[i][2] += fq*tizcoul; if (newton_pair || j < nlocal) { f[j][0] -= fx; f[j][1] -= fy; f[j][2] -= fz; torque[j][0] += fq*tjxcoul; torque[j][1] += fq*tjycoul; torque[j][2] += fq*tjzcoul; } if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { ecoul = (1.0-sqrt(rsq/cut_coulsq[itype][jtype])); ecoul *= ecoul; ecoul *= qtmp * q[j] * rinv; if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); if (mu[i][3] > 0.0 && q[j] != 0.0) ecoul += -q[j] * r3inv * pqfac * pidotr; if (mu[j][3] > 0.0 && qtmp != 0.0) ecoul += qtmp * r3inv * qpfac * pjdotr; ecoul *= factor_coul*qqrd2e*scale[itype][jtype]; } else ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype])+ rcutlj6inv*(6*lj3[itype][jtype]*rcutlj6inv-3*lj4[itype][jtype])* rsq*rcutlj2inv+ rcutlj6inv*(-7*lj3[itype][jtype]*rcutlj6inv+4*lj4[itype][jtype]); evdwl *= factor_lj; } else evdwl = 0.0; } if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, evdwl,ecoul,fx,fy,fz,delx,dely,delz); } } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::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(cut_lj,n+1,n+1,"pair:cut_lj"); memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); memory->create(cut_coul,n+1,n+1,"pair:cut_coul"); memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(lj1,n+1,n+1,"pair:lj1"); memory->create(lj2,n+1,n+1,"pair:lj2"); memory->create(lj3,n+1,n+1,"pair:lj3"); memory->create(lj4,n+1,n+1,"pair:lj4"); memory->create(scale,n+1,n+1,"pair:scale"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::settings(int narg, char **arg) { if (narg < 1 || narg > 2) error->all(FLERR,"Incorrect args in pair_style command"); if (strcmp(update->unit_style,"electron") == 0) error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles"); cut_lj_global = force->numeric(FLERR,arg[0]); if (narg == 1) cut_coul_global = cut_lj_global; else cut_coul_global = force->numeric(FLERR,arg[1]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) { cut_lj[i][j] = cut_lj_global; cut_coul[i][j] = cut_coul_global; } } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::coeff(int narg, char **arg) { - if (narg < 4 || narg > 7) + if (narg < 4 || narg > 8) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); - double scale_one = 1.0; - if (narg >= 5) scale_one = force->numeric(FLERR,arg[4]); - double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; - if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]); - if (narg == 7) cut_coul_one = force->numeric(FLERR,arg[6]); + double scale_one = 1.0; + int iarg = 4; + + if ((narg > iarg) && (strcmp(arg[iarg],"scale") != 0)) { + cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[iarg]); + ++iarg; + } + if ((narg > iarg) && (strcmp(arg[iarg],"scale") != 0)) { + cut_coul_one = force->numeric(FLERR,arg[iarg]); + ++iarg; + } + if (narg > iarg) { + if (strcmp(arg[iarg],"scale") == 0) { + scale_one = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else error->all(FLERR,"Incorrect args for pair coefficients"); + } + if (iarg != narg) + error->all(FLERR,"Incorrect args for pair coefficients"); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; cut_coul[i][j] = cut_coul_one; setflag[i][j] = 1; scale[i][j] = scale_one; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::init_style() { if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag) error->all(FLERR,"Pair dipole/sf requires atom attributes q, mu, torque"); neighbor->request(this,instance_me); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJSFDipoleSF::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]); } double cut = MAX(cut_lj[i][j],cut_coul[i][j]); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); cut_ljsq[j][i] = cut_ljsq[i][j]; cut_coulsq[j][i] = cut_coulsq[i][j]; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; scale[j][i] = scale[i][j]; return cut; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&epsilon[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut_lj[i][j],sizeof(double),1,fp); fwrite(&cut_coul[i][j],sizeof(double),1,fp); fwrite(&scale[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&epsilon[i][j],sizeof(double),1,fp); fread(&sigma[i][j],sizeof(double),1,fp); fread(&cut_lj[i][j],sizeof(double),1,fp); fread(&cut_coul[i][j],sizeof(double),1,fp); fread(&scale[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::write_restart_settings(FILE *fp) { fwrite(&cut_lj_global,sizeof(double),1,fp); fwrite(&cut_coul_global,sizeof(double),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJSFDipoleSF::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&cut_lj_global,sizeof(double),1,fp); fread(&cut_coul_global,sizeof(double),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } // PairLJSFDipoleSF: calculation of force is missing (to be implemented) double PairLJSFDipoleSF::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv; double pdotp,pidotr,pjdotr,pre1,delx,dely,delz; double rinv, r3inv,r5inv, rcutlj2inv, rcutcoul2inv,rcutlj6inv; double qtmp,xtmp,ytmp,ztmp,bfac,pqfac,qpfac, ecoul, evdwl; double **x = atom->x; double *q = atom->q; double **mu = atom->mu; + if (!warn_single) { + warn_single = 1; + if (comm->me == 0) { + error->warning(FLERR,"Single method for lj/sf/dipole/sf does not compute forces"); + } + } qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; r2inv = 1.0/rsq; rinv = sqrt(r2inv); fforce = 0.0; if (rsq < cut_coulsq[itype][jtype]) { delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; // if (qtmp != 0.0 && q[j] != 0.0) { // pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]); // forcecoulx += pre1*delx; // forcecouly += pre1*dely; // forcecoulz += pre1*delz; // } if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) + 3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv; } if (mu[i][3] > 0.0 && q[j] != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); } if (mu[j][3] > 0.0 && qtmp != 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; qpfac = 1.0 - 3.0*rsq*rcutcoul2inv + 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); } } if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; rcutlj2inv = 1.0 / cut_ljsq[itype][jtype]; rcutlj6inv = rcutlj2inv * rcutlj2inv * rcutlj2inv; } double eng = 0.0; if (rsq < cut_coulsq[itype][jtype]) { ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])); ecoul *= ecoul; ecoul *= qtmp * q[j] * rinv; if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); if (mu[i][3] > 0.0 && q[j] != 0.0) ecoul += -q[j] * r3inv * pqfac * pidotr; if (mu[j][3] > 0.0 && qtmp != 0.0) ecoul += qtmp * r3inv * qpfac * pjdotr; ecoul *= factor_coul*force->qqrd2e*scale[itype][jtype]; eng += ecoul; } if (rsq < cut_ljsq[itype][jtype]) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype])+ rcutlj6inv*(6*lj3[itype][jtype]*rcutlj6inv-3*lj4[itype][jtype])* rsq*rcutlj2inv+ rcutlj6inv*(-7*lj3[itype][jtype]*rcutlj6inv+4*lj4[itype][jtype]); eng += evdwl*factor_lj; } return eng; } /* ---------------------------------------------------------------------- */ void *PairLJSFDipoleSF::extract(const char *str, int &dim) { dim = 2; if (strcmp(str,"epsilon") == 0) return (void *) epsilon; if (strcmp(str,"sigma") == 0) return (void *) sigma; if (strcmp(str,"scale") == 0) return (void *) scale; return NULL; }