ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

New pair style lj/relres

Browse Source
This commit is contained in:
Mark Chaimovich
2021-02-09 07:29:38 -05:00
parent 222d842b45
commit 69f5d840df
10 changed files with 41836 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_pair.rst
View File

@ -163,6 +163,7 @@ OPT.
* :doc:`lj/long/dipole/long <pair_dipole>`
* :doc:`lj/long/tip4p/long (o) <pair_lj_long>`
* :doc:`lj/mdf <pair_mdf>`
* :doc:`lj/relres <pair_lj_relres>`
* :doc:`lj/sdk (gko) <pair_sdk>`
* :doc:`lj/sdk/coul/long (go) <pair_sdk>`
* :doc:`lj/sdk/coul/msm (o) <pair_sdk>`

140
doc/src/pair_lj_relres.rst Normal file
View File

@ -0,0 +1,140 @@
.. index:: pair_style lj/relres
pair_style lj/relres command
============================
Syntax
""""""
.. code-block:: LAMMPS
pair_style lj/relres Rsi Rso Rci Rco
* Rsi = inner switching distance - boundary up to which LJ potential of fine-grained model is applied (distance units)
* Rso = outer switching distance - boundary beyond which LJ potential of coarse-grained model is applied (distance units)
* Rci = inner cutoff beyond which force smoothing is applied (distance units)
* Rco = outer cutoff for lj/relres interactions (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style lj/relres 4.0 5.0 8.0 10.0
pair_coeff * * 0.5 1.0 0.0 1.0
pair_coeff 1 1 1.5 1.2 3.0 1.2 3.0 3.5 6.0 7.0
Description
"""""""""""
Style *lj/relres* computes a LJ interaction with RelRes methodology developed by :ref:`Chaimovich at al.<Chaimovich1>`
This methodology applies fine-grained model between near neighbors (up to :math:`r_{si}` boundary) and a simplified coarse-grained model
for far neighbors (beyond :math:`r_{so}` boundary) allowing significant improvement in computational efficiency while preserving correctness
of simulation results.
.. math::
E & = & \left\{\begin{array}{lr}
4 \epsilon^{FG} \left[ \left(\frac{\sigma^{FG}}{r}\right)^{12} - \left(\frac{\sigma^{FG}}{r}\right)^6 \right]-G_{si}, & r< r_{si} \\
\sum_{m=0}^{4} C_{sm}\left(r-r_{si}\right)^m-G_{so} , & r_{si}\leq r< r_{so} \\
4 \epsilon^{CG} \left[ \left(\frac{\sigma^{CG}}{r}\right)^{12} - \left(\frac{\sigma^{CG}}{r}\right)^6 \right]-G_c, & r_{so}\leq r<r_{ci} \\
\sum_{m=0}^{4} C_{cm}\left(r-r_{ci}\right)^m -G_c, & r_{ci}\leq r< r_{co} \\
0, & r\geq r_{co}\end{array}\right. \\
& &
Between :math:`r_{si}` and :math:`r_{so}` the polynomial smoothing is applied in a way that the force and its 1st derivative are not discontinued
at switching between fine- and coarse-grained potentials (between :math:`r_{si}` and :math:`r_{so}`) and at cutoff (between :math:`r_{ci}` and :math:`r_{co}`).
The corresponding polynomial coefficients :math:`C_{sm}` and :math:`C_{cm}` and shifting constants :math:`G_{si}`, :math:`G_{so}` and :math:`G_{c}` are computed by LAMMPS accordingly.
To avoid smoothing, the inner switching distance :math:`r_{si}` parameter should be set equal to the outer switching distance :math:`r_{so}` parameter
(:math:`r_{si}=r_{so}`). Similarly, to avoid smoothing at cutoff, inner and outer cutoff parameters should be set equal (:math:`r_{ci}=r_{co}`).
Details can be found in :ref:`(Chaimovich) <Chaimovich2>`.
.. note::
Energy and force resulting from this methodology can be plotted via the
:doc:`pair_write <pair_write>` command to see the effect.
In implementation of *lj/relres* style, atoms are grouped in the way that one of the atoms in the group plays the role of a coarse-grained site for the calculation
of interactions beyond :math:`r_{so}` distance while continuing to play the role of a fine-grained site for shorter distances.
This atom must be defined as a different atom type. Other atoms in the group participate in the fine-grained interactions only.
The following coefficients must be defined for each pair of atom
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands, or by mixing as will be described below:
* :math:`\epsilon^{FG}` (energy units)
* :math:`\sigma^{FG}` (distance units)
* :math:`\epsilon^{CG}` (energy units)
* :math:`\sigma^{CG}` (distance units)
For atom types that are used as fine-grained sites only, :math:`\epsilon^{CG}` must be set to 0 (zero).
For atom types that are used as coarse-grained sites only (if any), :math:`\epsilon^{FG}` must be set to 0 (zero).
Additional parameters can be defined to specify different :math:`r_{si}`, :math:`r_{so}`, :math:`r_{ci}`, :math:`r_{co}` for a particular set of atom types:
* :math:`r_{si}` (distance units)
* :math:`r_{so}` (distance units)
* :math:`r_{ci}` (distance units)
* :math:`r_{co}` (distance units)
These parameters are optional and they are used to override global values defined in the pair_style command.
If this override option is used, all four values must be specified. If not specified, the global values for :math:`r_{si}`, :math:`r_{so}`, :math:`r_{ci}`, and :math:`r_{co}` are used.
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J with I != J, the :math:`\epsilon^{FG}`, :math:`\sigma^{FG}`, :math:`\epsilon^{CG}`, :math:`\sigma^{CG}`, :math:`r_{si}`, :math:`r_{so}`, :math:`r_{ci}`, and :math:`r_{co}`
parameters for this pair style can be mixed, if not defined explicitly.
All parameters are mixed according to the pair_modify mix option. The
default mix value is *geometric*\ , and it is recommended to use with this *lj/relres* style.
See the "pair_modify" command for details.
This pair style supports the :doc:`pair_modify <pair_modify>` shift
option for the energy of the pair interaction. It is recommended to set this option to *yes*\ .
Otherwise, the shifting constant :math:`G_{c}` is set to zero. Constants :math:`G_{si}` and :math:`G_{so}` are not impacted by this option.
The :doc:`pair_modify <pair_modify>` table option is not relevant
for this pair style.
This pair style does not support the :doc:`pair_modify <pair_modify>`
tail option for adding long-range tail corrections to energy and
pressure, since the energy of the pair interaction is smoothed to 0.0
at the cutoff.
This pair style writes its information to :doc:`binary restart files <restart>`, 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
:doc:`run_style respa <run_style>` command. It does not support the
*inner*\ , *middle*\ , *outer* keywords.
----------
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`
Default
"""""""
none
----------
.. _Chaimovich1:
**(Chaimovich at al.)** A.Chaimovich, C. Peter and K. Kremer, J. Chem. Phys. 143, 243107
(2015).
.. _Chaimovich2:
**(Chaimovich)** M.Chaimovich and A. Chaimovich, J. Chem. Theory Comput. 17, 1045-1059
(2021).

1
doc/src/pair_style.rst
View File

@ -227,6 +227,7 @@ accelerated styles exist.
* :doc:`lj/long/dipole/long <pair_dipole>` - long-range LJ and long-range point dipoles
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water
* :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function
* :doc:`lj/relres <pair_lj_relres>` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) <Chaimovich2>`.
* :doc:`lj/sdk <pair_sdk>` - LJ for SDK coarse-graining
* :doc:`lj/sdk/coul/long <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulomb
* :doc:`lj/sdk/coul/msm <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulomb via MSM

37075
examples/relres/Data.22DMH.in.relres Normal file
View File

File diff suppressed because it is too large Load Diff

5
examples/relres/README Normal file
View File

@ -0,0 +1,5 @@
The input script in.22DMH.relres provide example of simulation using Relative Resolution (RelRes) potential.
In this example 2,2-Dimethylhexane is selected as simulated substance to give complete view of the RelRes utilization.
This script uses data file Data.22DMH.in.relres consisting 8000 molecule of the substance. It also generates RelRes potential for selected atom types.

31
examples/relres/in.22DMH.relres Normal file
View File

@ -0,0 +1,31 @@
units si
atom_style molecular
boundary p p p
dielectric 1
special_bonds lj/coul 0.0 0.0 0.5
pair_style lj/relres 0.675e-9 .725e-9 1.2e-9 1.4e-9
bond_style harmonic
angle_style harmonic
dihedral_style fourier
read_data Data.22DMH.in.relres
pair_coeff 6 6 1.21585e-21 0.3905e-9 0 0.3905e-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 4 4 0.819828e-21 0.3905e-9 0 0.3905e-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 2 2 0.819828e-21 0.3905e-9 8.48872E-21 0.3905E-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 5 5 1.00742E-21 0.396E-9 0 0.396E-9
pair_coeff 3 3 0.819828e-21 0.3905e-9 0 0.3905e-9
pair_coeff 1 1 0.347385E-21 0.38E-9 20.2372E-21 0.39309E-9
pair_modify shift yes
neighbor 2.0e-10 bin
neigh_modify every 2 delay 4 check yes
timestep 1.0e-15
thermo 100
thermo_style custom step temp press pe ke etotal epair emol vol
fix 2 all nvt temp 290 290 2.0e-13
run 1000
#write_data Data.22DMH.out.relres pair ij
pair_write 1 1 1201 r 0.2e-9 1.4e-9 potential.relres LJ11
pair_write 2 2 1201 r 0.2e-9 1.4e-9 potential.relres LJ22
pair_write 1 2 1201 r 0.2e-9 1.4e-9 potential.relres LJ12

121
examples/relres/log.relres Normal file
View File

@ -0,0 +1,121 @@
LAMMPS (24 Dec 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
using 1 OpenMP thread(s) per MPI task
units si
atom_style molecular
boundary p p p
dielectric 1
special_bonds lj/coul 0.0 0.0 0.5
pair_style lj/relres 0.675e-9 .725e-9 1.2e-9 1.4e-9
bond_style harmonic
angle_style harmonic
dihedral_style fourier
read_data Data.22DMH.in.relres
Reading data file ...
orthogonal box = (3.7421629e-10 3.7421629e-10 3.7421629e-10) to (6.8257837e-09 6.8257837e-09 6.8257837e-09)
1 by 1 by 2 MPI processor grid
reading atoms ...
8000 atoms
reading velocities ...
8000 velocities
scanning bonds ...
4 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
3 = max dihedrals/atom
reading bonds ...
7000 bonds
reading angles ...
9000 angles
reading dihedrals ...
5000 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0.5
special bond factors coul: 0 0 0.5
4 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
5 = max # of 1-4 neighbors
7 = max # of special neighbors
special bonds CPU = 0.003 seconds
read_data CPU = 0.062 seconds
pair_coeff 6 6 1.21585e-21 0.3905e-9 0 0.3905e-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 4 4 0.819828e-21 0.3905e-9 0 0.3905e-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 2 2 0.819828e-21 0.3905e-9 8.48872E-21 0.3905E-9 0.575e-9 0.625e-9 1.2e-9 1.4e-9
pair_coeff 5 5 1.00742E-21 0.396E-9 0 0.396E-9
pair_coeff 3 3 0.819828e-21 0.3905e-9 0 0.3905e-9
pair_coeff 1 1 0.347385E-21 0.38E-9 20.2372E-21 0.39309E-9
pair_modify shift yes
neighbor 2.0e-10 bin
neigh_modify every 2 delay 4 check yes
timestep 1.0e-15
thermo 100
thermo_style custom step temp press pe ke etotal epair emol vol
fix 2 all nvt temp 290 290 2.0e-13
run 1000
Neighbor list info ...
update every 2 steps, delay 4 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.6e-09
ghost atom cutoff = 1.6e-09
binsize = 8e-10, bins = 9 9 9
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/relres, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 19.53 | 19.53 | 19.53 Mbytes
Step Temp Press PotEng KinEng TotEng E_pair E_mol Volume
0 286.85659 -21390710 2.6345656e-17 4.7519899e-17 7.3865555e-17 -3.3260066e-17 5.9605722e-17 2.685318e-25
100 292.25165 10716172 2.5245584e-17 4.841363e-17 7.3659214e-17 -3.2561964e-17 5.7807548e-17 2.685318e-25
200 291.6011 48774461 2.4897987e-17 4.8305863e-17 7.320385e-17 -3.3268705e-17 5.8166692e-17 2.685318e-25
300 291.50656 37655969 2.4389062e-17 4.8290201e-17 7.2679262e-17 -3.3428236e-17 5.7817297e-17 2.685318e-25
400 287.23427 25920755 2.4747225e-17 4.7582464e-17 7.2329689e-17 -3.3065908e-17 5.7813133e-17 2.685318e-25
500 288.56911 -9297451 2.4379025e-17 4.7803591e-17 7.2182615e-17 -3.3515426e-17 5.7894451e-17 2.685318e-25
600 291.82949 20083719 2.3686904e-17 4.8343696e-17 7.2030599e-17 -3.3468666e-17 5.7155569e-17 2.685318e-25
700 290.64445 60535932 2.3704156e-17 4.8147386e-17 7.1851542e-17 -3.3299994e-17 5.700415e-17 2.685318e-25
800 293.01243 38119194 2.3163674e-17 4.8539659e-17 7.1703333e-17 -3.3641284e-17 5.6804958e-17 2.685318e-25
900 289.1191 32514067 2.3608264e-17 4.7894701e-17 7.1502965e-17 -3.3787865e-17 5.7396129e-17 2.685318e-25
1000 292.45834 -714652.68 2.2885394e-17 4.8447871e-17 7.1333265e-17 -3.3864138e-17 5.6749532e-17 2.685318e-25
Loop time of 7.07676 on 2 procs for 1000 steps with 8000 atoms
99.7% CPU use with 2 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.1275 | 5.1852 | 5.243 | 2.5 | 73.27
Bond | 1.0056 | 1.0072 | 1.0087 | 0.2 | 14.23
Neigh | 0.38261 | 0.38263 | 0.38266 | 0.0 | 5.41
Comm | 0.25513 | 0.31266 | 0.37018 | 10.3 | 4.42
Output | 0.00042391 | 0.00044036 | 0.00045681 | 0.0 | 0.01
Modify | 0.14501 | 0.14643 | 0.14785 | 0.4 | 2.07
Other | | 0.04219 | | | 0.60
Nlocal: 4000.00 ave 4000 max 4000 min
Histogram: 2 0 0 0 0 0 0 0 0 0
Nghost: 13887.5 ave 13896 max 13879 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Neighs: 217142.0 ave 217471 max 216812 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Total # of neighbors = 434283
Ave neighs/atom = 54.285375
Ave special neighs/atom = 5.2500000
Neighbor list builds = 13
Dangerous builds = 0
#write_data Data.22DMH.out.relres pair ij
pair_write 1 1 1201 r 0.2e-9 1.4e-9 potential.relres LJ11
Creating table file potential.relres with DATE: 2021-02-05
pair_write 2 2 1201 r 0.2e-9 1.4e-9 potential.relres LJ22
Appending to table file potential.relres with DATE: 2021-02-05
pair_write 1 2 1201 r 0.2e-9 1.4e-9 potential.relres LJ12
Appending to table file potential.relres with DATE: 2021-02-05
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:07

3619
examples/relres/potential.relres Normal file
View File

File diff suppressed because it is too large Load Diff

769
src/pair_lj_relres.cpp Normal file
View File

@ -0,0 +1,769 @@
/* ----------------------------------------------------------------------
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: Mark Chaimovich(RSM) mark.chaimovich@russianschool.com
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include "pair_lj_relres.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
using namespace LAMMPS_NS;
static const char cite_relres[] =
"RelRes:\n\n"
"@Article{Chaimovich1,\n"
" author = {A. Chaimovich, C. Peter, K. Kremer},\n"
" title = {Relative resolution: A hybrid formalism for fluid mixtures},\n"
" journal = {J.~Chem.~Phys.},\n"
" year = 2015,\n"
" volume = 143,\n"
" pages = {243107}\n"
"@Article{Chaimovich2,\n"
" author = {M. Chaimovich, A. Chaimovich},\n"
" title = {Relative Resolution: A Computationally Efficient Implementation in LAMMPS},\n"
" journal = {J.~Chem.~Theory~Comput.},\n"
" year = 2021,\n"
" volume = 17,\n"
" pages = {1045--1059}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairLJRelRes::PairLJRelRes(LAMMPS *lmp) : Pair(lmp)
{
if (lmp->citeme) lmp->citeme->add(cite_relres);
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairLJRelRes::~PairLJRelRes()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cutfsq);
memory->destroy(cut);
memory->destroy(cut_inner);
memory->destroy(cut_inner_sq);
memory->destroy(cutf);
memory->destroy(cutf_inner);
memory->destroy(cutf_inner_sq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(epsilonf);
memory->destroy(sigmaf);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(ljsw0);
memory->destroy(ljsw1);
memory->destroy(ljsw2);
memory->destroy(ljsw3);
memory->destroy(ljsw4);
memory->destroy(ljf1);
memory->destroy(ljf2);
memory->destroy(ljf3);
memory->destroy(ljf4);
memory->destroy(ljswf0);
memory->destroy(ljswf1);
memory->destroy(ljswf2);
memory->destroy(ljswf3);
memory->destroy(ljswf4);
memory->destroy(offset);
memory->destroy(offsetsm);
memory->destroy(offsetsp);
}
}
/* ---------------------------------------------------------------------- */
void PairLJRelRes::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj;
double r,t,tsq,fskin;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
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];
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)];
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;
if (rsq < cutf_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(ljf1[itype][jtype]*r6inv-ljf2[itype][jtype]);
}
else
if (rsq < cutfsq[itype][jtype]) {
r = sqrt(rsq);
t = r - cutf_inner[itype][jtype];
tsq = t*t;
fskin = ljswf1[itype][jtype]+ljswf2[itype][jtype]*t+
ljswf3[itype][jtype]*tsq+ljswf4[itype][jtype]*tsq*t;
forcelj = fskin*r;
}
else
if (rsq < cut_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
}
else {
r = sqrt(rsq);
t = r-cut_inner[itype][jtype];
tsq = t*t;
fskin = ljsw1[itype][jtype]+ljsw2[itype][jtype]*t+
ljsw3[itype][jtype]*tsq+ljsw4[itype][jtype]*tsq*t;
forcelj = fskin*r;
}
fpair = factor_lj*forcelj*r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
if (rsq < cutf_inner_sq[itype][jtype])
evdwl = r6inv*(ljf3[itype][jtype]*r6inv-
ljf4[itype][jtype])-offsetsm[itype][jtype];
else
if (rsq < cutfsq[itype][jtype])
evdwl = ljswf0[itype][jtype]-ljswf1[itype][jtype]*t-
ljswf2[itype][jtype]*tsq/2.0-ljswf3[itype][jtype]*tsq*t/3.0-
ljswf4[itype][jtype]*tsq*tsq/4.0-offsetsp[itype][jtype];
else
if (rsq < cut_inner_sq[itype][jtype])
evdwl = r6inv*(lj3[itype][jtype]*r6inv-
lj4[itype][jtype])-offset[itype][jtype];
else
evdwl = ljsw0[itype][jtype]-ljsw1[itype][jtype]*t-
ljsw2[itype][jtype]*tsq/2.0-ljsw3[itype][jtype]*tsq*t/3.0-
ljsw4[itype][jtype]*tsq*tsq/4.0-offset[itype][jtype];
evdwl *= factor_lj;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJRelRes::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(cutfsq,n+1,n+1,"pair:cutfsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cutf,n+1,n+1,"pair:cutf");
memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
memory->create(cutf_inner,n+1,n+1,"pair:cutf_inner");
memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq");
memory->create(cutf_inner_sq,n+1,n+1,"pair:cutf_inner_sq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(epsilonf,n+1,n+1,"pair:epsilonf");
memory->create(sigmaf,n+1,n+1,"pair:sigmaf");
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(ljf1,n+1,n+1,"pair:ljf1");
memory->create(ljf2,n+1,n+1,"pair:ljf2");
memory->create(ljf3,n+1,n+1,"pair:ljf3");
memory->create(ljf4,n+1,n+1,"pair:ljf4");
memory->create(ljsw0,n+1,n+1,"pair:ljsw0");
memory->create(ljsw1,n+1,n+1,"pair:ljsw1");
memory->create(ljsw2,n+1,n+1,"pair:ljsw2");
memory->create(ljsw3,n+1,n+1,"pair:ljsw3");
memory->create(ljsw4,n+1,n+1,"pair:ljsw4");
memory->create(ljswf0,n+1,n+1,"pair:ljswf0");
memory->create(ljswf1,n+1,n+1,"pair:ljswf1");
memory->create(ljswf2,n+1,n+1,"pair:ljswf2");
memory->create(ljswf3,n+1,n+1,"pair:ljswf3");
memory->create(ljswf4,n+1,n+1,"pair:ljswf4");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(offsetsp,n+1,n+1,"pair:offsetsp");
memory->create(offsetsm,n+1,n+1,"pair:offsetsm");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJRelRes::settings(int narg, char **arg)
{
if (narg != 4) error->all(FLERR,"Illegal pair_style command");
cutf_inner_global = utils::numeric(FLERR,arg[0],false,lmp);
cutf_global = utils::numeric(FLERR,arg[1],false,lmp);
cut_inner_global = utils::numeric(FLERR,arg[2],false,lmp);
cut_global = utils::numeric(FLERR,arg[3],false,lmp);
if (cut_inner_global <= 0.0 || cut_inner_global > cut_global)
error->all(FLERR,"Illegal pair_style command");
if (cutf_inner_global <= 0.0 || cutf_inner_global > cutf_global)
error->all(FLERR,"Illegal pair_style command");
if (cutf_global > cut_inner_global)
error->all(FLERR,"Illegal pair_style command");
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_inner[i][j] = cut_inner_global;
cut[i][j] = cut_global;
cutf_inner[i][j] = cutf_inner_global;
cutf[i][j] = cutf_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJRelRes::coeff(int narg, char **arg)
{
if (narg != 6 && narg != 10)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
double epsilonf_one = utils::numeric(FLERR,arg[2],false,lmp);
double sigmaf_one = utils::numeric(FLERR,arg[3],false,lmp);
double epsilon_one = utils::numeric(FLERR,arg[4],false,lmp);
double sigma_one = utils::numeric(FLERR,arg[5],false,lmp);
double cut_inner_one = cut_inner_global;
double cut_one = cut_global;
double cutf_inner_one = cutf_inner_global;
double cutf_one = cutf_global;
if (narg == 10) {
cutf_inner_one = utils::numeric(FLERR,arg[6],false,lmp);
cutf_one = utils::numeric(FLERR,arg[7],false,lmp);
cut_inner_one = utils::numeric(FLERR,arg[8],false,lmp);
cut_one = utils::numeric(FLERR,arg[9],false,lmp);
}
if (cut_inner_one <= 0.0 || cut_inner_one > cut_one)
error->all(FLERR,"Incorrect args for pair coefficients");
if (cutf_inner_one <= 0.0 || cutf_inner_one > cutf_one)
error->all(FLERR,"Incorrect args for pair coefficients");
if (cutf_one > cut_inner_one)
error->all(FLERR,"Incorrect args for pair coefficients");
if (epsilon_one == 0.0) { //set cutoff for fg interactions
cut_inner_one = cutf_one;
cut_one = cutf_one;
}
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;
epsilonf[i][j] = epsilonf_one;
sigmaf[i][j] = sigmaf_one;
cut_inner[i][j] = cut_inner_one;
cut[i][j] = cut_one;
cutf_inner[i][j] = cutf_inner_one;
cutf[i][j] = cutf_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJRelRes::init_one(int i, int j)
{ double ljswc0,ljswc3,ljswc4;
// mixing rules:
// fg and cg - no mixing;
// fg and fg or fg anf hybrid - mixing fg parameters only
// cg and cg of cg and hybrid - mixing cg parameters only
// hybrid and hybrid - mixing fg and cg parameters
if (setflag[i][j] == 0) {
if (epsilon[i][i] == 0.0 && epsilonf[j][j] == 0.0 ||
epsilonf[i][i] == 0.0 && epsilon[j][j] == 0.0) { //no mixing
epsilon[i][j] = 0.0;
epsilonf[i][j] = 0.0;
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
sigmaf[i][j] = mix_distance(sigmaf[i][i],sigmaf[j][j]);
cut_inner[i][j] = cutf[i][j] = cutf_inner[i][j] = cut[i][j] = 0.0;
}
else
if (epsilon[i][i] == 0.0 || epsilon[j][j] == 0.0) { // fg only
epsilon[i][j] = 0.0;
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
epsilonf[i][j] = mix_energy(epsilonf[i][i],epsilonf[j][j],
sigmaf[i][i],sigmaf[j][j]);
sigmaf[i][j] = mix_distance(sigmaf[i][i],sigmaf[j][j]);
cutf_inner[i][j] = mix_distance(cutf_inner[i][i],cutf_inner[j][j]);
cutf[i][j] = mix_distance(cutf[i][i],cutf[j][j]);
cut_inner[i][j] = cutf[i][j];
cut[i][j] = cutf[i][j];
}
else
if (epsilonf[i][i] == 0.0 || epsilonf[j][j] == 0.0) { // cg only
epsilonf[i][j] = 0.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]);
sigmaf[i][j] = mix_distance(sigmaf[i][i],sigmaf[j][j]);
cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
cutf_inner[i][j] = mix_distance(cutf_inner[i][i],cutf_inner[j][j]);
cutf[i][j] = mix_distance(cutf[i][i],cutf[j][j]);
}
else { // fg and cg
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]);
epsilonf[i][j] = mix_energy(epsilonf[i][i],epsilonf[j][j],
sigmaf[i][i],sigmaf[j][j]);
sigmaf[i][j] = mix_distance(sigmaf[i][i],sigmaf[j][j]);
cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
cutf_inner[i][j] = mix_distance(cutf_inner[i][i],cutf_inner[j][j]);
cutf[i][j] = mix_distance(cutf[i][i],cutf[j][j]);
}
}
cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j];
cutf_inner_sq[i][j] = cutf_inner[i][j]*cutf_inner[i][j];
cutfsq[i][j] = cutf[i][j]*cutf[i][j];
if (epsilon[i][j] != 0) { // cg or fg+cg (cut coefficients)
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);
if (cut_inner[i][j] != cut[i][j]) {
double r6inv = 1.0/pow(cut_inner[i][j],6.0);
double t = cut[i][j] - cut_inner[i][j];
double tsq = t*t;
double ratio = sigma[i][j] / cut_inner[i][j];
ljsw0[i][j] = 4.0*epsilon[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
ljsw1[i][j] = r6inv*(lj1[i][j]*r6inv-lj2[i][j]) / cut_inner[i][j];
ljsw2[i][j] = -r6inv * (13.0*lj1[i][j]*r6inv - 7.0*lj2[i][j]) /
cut_inner_sq[i][j];
ljsw3[i][j] = -(3.0/tsq) * (ljsw1[i][j] + 2.0/3.0*ljsw2[i][j]*t);
ljsw4[i][j] = -1.0/(3.0*tsq) * (ljsw2[i][j] + 2.0*ljsw3[i][j]*t);
if (offset_flag)
offset[i][j] = ljsw0[i][j] - ljsw1[i][j]*t - ljsw2[i][j]*tsq/2.0 -
ljsw3[i][j]*tsq*t/3.0 - ljsw4[i][j]*tsq*tsq/4.0;
else offset[i][j] = 0.0;
}
else {
ljsw0[i][j] = 0.0;
ljsw1[i][j] = 0.0;
ljsw2[i][j] = 0.0;
ljsw3[i][j] = 0.0;
ljsw4[i][j] = 0.0;
double ratio = sigma[i][j] / cut_inner[i][j];
if (offset_flag)
offset[i][j] = 4.0*epsilon[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
else offset[i][j] = 0.0;
}
}
else {
ljsw0[i][j] = 0.0;
ljsw1[i][j] = 0.0;
ljsw2[i][j] = 0.0;
ljsw3[i][j] = 0.0;
ljsw4[i][j] = 0.0;
lj1[i][j] = 0.0;
lj2[i][j] = 0.0;
lj3[i][j] = 0.0;
lj4[i][j] = 0.0;
offset[i][j] = 0.0;
}
if (epsilonf[i][j] != 0 ) { // fg (cut=cutf coefficients)
ljf1[i][j] = 48.0 * epsilonf[i][j] * pow(sigmaf[i][j],12.0);
ljf2[i][j] = 24.0 * epsilonf[i][j] * pow(sigmaf[i][j],6.0);
ljf3[i][j] = 4.0 * epsilonf[i][j] * pow(sigmaf[i][j],12.0);
ljf4[i][j] = 4.0 * epsilonf[i][j] * pow(sigmaf[i][j],6.0);
if (cutf_inner[i][j] != cutf[i][j]) {
double r6inv = 1.0/pow(cutf_inner[i][j],6.0);
double t = cutf[i][j] - cutf_inner[i][j];
double tsq = t*t;
double ratio = sigmaf[i][j] / cutf_inner[i][j];
ljswf0[i][j] = 4.0*epsilonf[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
ljswf1[i][j] = r6inv*(ljf1[i][j]*r6inv-ljf2[i][j]) / cutf_inner[i][j];
ljswf2[i][j] = -r6inv * (13.0*ljf1[i][j]*r6inv - 7.0*ljf2[i][j]) /
cutf_inner_sq[i][j];
ljswf3[i][j] = -(3.0/tsq) * (ljswf1[i][j] + 2.0/3.0*ljswf2[i][j]*t);
ljswf4[i][j] = -1.0/(3.0*tsq) * (ljswf2[i][j] + 2.0*ljswf3[i][j]*t);
offsetsp[i][j] = ljswf0[i][j] - ljswf1[i][j]*t - ljswf2[i][j]*tsq/2.0-
ljswf3[i][j]*tsq*t/3.0 - ljswf4[i][j]*tsq*tsq/4.0;
}
else {
ljswf0[i][j] = 0.0;
ljswf1[i][j] = 0.0;
ljswf2[i][j] = 0.0;
ljswf3[i][j] = 0.0;
ljswf4[i][j] = 0.0;
double ratio = sigmaf[i][j] / cutf_inner[i][j];
offsetsp[i][j] = 4.0*epsilonf[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
}
}
else {
ljswf0[i][j] = 0.0;
ljswf1[i][j] = 0.0;
ljswf2[i][j] = 0.0;
ljswf3[i][j] = 0.0;
ljswf4[i][j] = 0.0;
ljf4[i][j] = 0.0;
ljf1[i][j] = 0.0;
ljf2[i][j] = 0.0;
ljf3[i][j] = 0.0;
offsetsp[i][j] = 0.0;
}
if (epsilon[i][j] != 0) { // cg or fg+cg (cutf coefficients)
if (cutf_inner[i][j] != cutf[i][j]) {
double r2inv = 1.0/pow(cutf[i][j],2.0);
double r6inv = r2inv * r2inv * r2inv;
double t = cutf[i][j] - cutf_inner[i][j];
double tsq = t*t;
double tsqinv = 1.0/tsq;
double ratio = sigma[i][j] / cutf[i][j];
double Et = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
double Ft = r6inv * (lj1[i][j] * r6inv - lj2[i][j]) / cutf[i][j];
double dFt = -r6inv * (13.0*lj1[i][j]*r6inv - 7.0*lj2[i][j]) * r2inv;
double A = Ft + dFt * t / 3.0;
ljswc3 = 3.0 * A * tsqinv;
ljswc4 = -(2.0 * Ft + dFt * t) * tsqinv / t;
ljswc0 = Et + ljswc3 * t * tsq /3.0 + ljswc4 * tsq * tsq / 4.0;
offsetsm[i][j] = ljswc0;
}
else {
ljswc0 = 0.0;
ljswc3 = 0.0;
ljswc4 = 0.0;
double ratio = sigma[i][j] / cutf_inner[i][j];
offsetsm[i][j] = 4.0*epsilon[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
}
}
else {
ljswc0 = 0.0;
ljswc3 = 0.0;
ljswc4 = 0.0;
offsetsm[i][j] = 0.0;
}
// combine cutf coefficients
ljswf0[i][j] += ljswc0;
ljswf3[i][j] += ljswc3;
ljswf4[i][j] += ljswc4;
// combine shifting constants
offsetsp[i][j] += offset[i][j];
offsetsm[i][j] = offsetsp[i][j] - offsetsm[i][j];
if (i !=j) {
cut[j][i] = cut[i][j];
cutsq[j][i] = cutsq[i][j];
cutf[j][i] = cutf[i][j];
cutfsq[j][i] = cutfsq[i][j];
cut_inner[j][i] = cut_inner[i][j];
cut_inner_sq[j][i] = cut_inner_sq[i][j];
cutf_inner[j][i] = cutf_inner[i][j];
cutf_inner_sq[j][i] = cutf_inner_sq[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];
ljsw0[j][i] = ljsw0[i][j];
ljsw1[j][i] = ljsw1[i][j];
ljsw2[j][i] = ljsw2[i][j];
ljsw3[j][i] = ljsw3[i][j];
ljsw4[j][i] = ljsw4[i][j];
offset[j][i] = offset[i][j];
ljf1[j][i] = ljf1[i][j];
ljf2[j][i] = ljf2[i][j];
ljf3[j][i] = ljf3[i][j];
ljf4[j][i] = ljf4[i][j];
ljswf0[j][i] = ljswf0[i][j];
ljswf1[j][i] = ljswf1[i][j];
ljswf2[j][i] = ljswf2[i][j];
ljswf3[j][i] = ljswf3[i][j];
ljswf4[j][i] = ljswf4[i][j];
offsetsp[j][i] = offsetsp[i][j];
offsetsm[j][i] = offsetsm[i][j];
}
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJRelRes::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(&epsilonf[i][j],sizeof(double),1,fp);
fwrite(&sigmaf[i][j],sizeof(double),1,fp);
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cutf_inner[i][j],sizeof(double),1,fp);
fwrite(&cutf[i][j],sizeof(double),1,fp);
fwrite(&cut_inner[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJRelRes::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(&epsilonf[i][j],sizeof(double),1,fp);
fread(&sigmaf[i][j],sizeof(double),1,fp);
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cutf_inner[i][j],sizeof(double),1,fp);
fread(&cutf[i][j],sizeof(double),1,fp);
fread(&cut_inner[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilonf[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigmaf[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutf_inner[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutf[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJRelRes::write_restart_settings(FILE *fp)
{
fwrite(&cutf_inner_global,sizeof(double),1,fp);
fwrite(&cutf_global,sizeof(double),1,fp);
fwrite(&cut_inner_global,sizeof(double),1,fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJRelRes::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cutf_inner_global,sizeof(double),1,fp);
fread(&cutf_global,sizeof(double),1,fp);
fread(&cut_inner_global,sizeof(double),1,fp);
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cutf_inner_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutf_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairLJRelRes::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g %g %g %g\n",i,epsilonf[i][i],sigmaf[i][i],
epsilon[i][i],sigma[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairLJRelRes::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g %g %g %g %g %g\n",i,j,
epsilonf[i][j],sigmaf[i][j],epsilon[i][j],sigma[i][j],
cutf_inner[i][j],cutf[i][j],cut_inner[i][j],cut[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairLJRelRes::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r6inv,forcelj,philj,r,t,tsq,fskin;
r2inv = 1.0/rsq;
if (rsq < cutf_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(ljf1[itype][jtype]*r6inv-ljf2[itype][jtype]);
}
else
if (rsq < cutfsq[itype][jtype]) {
r = sqrt(rsq);
t = r - cutf_inner[itype][jtype];
tsq = t*t;
fskin = ljswf1[itype][jtype]+ljswf2[itype][jtype]*t+
ljswf3[itype][jtype]*tsq+ljswf4[itype][jtype]*tsq*t;
forcelj = fskin*r;
}
else
if (rsq < cut_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
}
else {
r = sqrt(rsq);
t = r - cut_inner[itype][jtype];
tsq = t*t;
fskin = ljsw1[itype][jtype] + ljsw2[itype][jtype]*t +
ljsw3[itype][jtype]*tsq + ljsw4[itype][jtype]*tsq*t;
forcelj = fskin*r;
}
fforce = factor_lj*forcelj*r2inv;
if (rsq < cutf_inner_sq[itype][jtype])
philj = r6inv*(ljf3[itype][jtype]*r6inv-
ljf4[itype][jtype])-offsetsm[itype][jtype];
else
if (rsq < cutfsq[itype][jtype])
philj = ljswf0[itype][jtype]-ljswf1[itype][jtype]*t-
ljswf2[itype][jtype]*tsq/2.0-ljswf3[itype][jtype]*tsq*t/3.0-
ljswf4[itype][jtype]*tsq*tsq/4.0-offsetsp[itype][jtype];
else
if (rsq < cut_inner_sq[itype][jtype])
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]) -
offset[itype][jtype];
else
philj = ljsw0[itype][jtype] - ljsw1[itype][jtype]*t -
ljsw2[itype][jtype]*tsq/2.0 - ljsw3[itype][jtype]*tsq*t/3.0 -
ljsw4[itype][jtype]*tsq*tsq/4.0 - offset[itype][jtype];
return factor_lj*philj;
}

74
src/pair_lj_relres.h Normal file
View File

@ -0,0 +1,74 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/relres,PairLJRelRes)
#else
#ifndef LMP_PAIR_LJ_RELRES_H
#define LMP_PAIR_LJ_RELRES_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJRelRes : public Pair {
public:
PairLJRelRes(class LAMMPS *);
virtual ~PairLJRelRes();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
double single(int, int, int, int, double, double, double, double &);
protected:
double cut_inner_global, cut_global, cutf_inner_global, cutf_global;
double **cut,**cut_inner,**cut_inner_sq,**cutf,**cutfsq,**cutf_inner,**cutf_inner_sq;
double **epsilon,**sigma;
double **epsilonf,**sigmaf;
double **lj1,**lj2,**lj3,**lj4;
double **ljf1,**ljf2,**ljf3,**ljf4;
double **ljsw0,**ljsw1,**ljsw2,**ljsw3,**ljsw4;
double **ljswf0,**ljswf1,**ljswf2,**ljswf3,**ljswf4;
double **offset,**offsetsp,**offsetsm;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/