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

merge with Axel changes

Browse Source
This commit is contained in:
Steve Plimpton
2022-04-19 09:14:56 -06:00
parent bae7fd74b0 c015851fbf
commit 44a4f3cff8
49 changed files with 1633 additions and 1594 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
cmake/CMakeLists.txt
View File

@ -167,6 +167,7 @@ option(CMAKE_VERBOSE_MAKEFILE "Generate verbose Makefiles" OFF)
set(STANDARD_PACKAGES
ADIOS
AMOEBA
ASPHERE
ATC
AWPMD

1
cmake/presets/all_off.cmake
View File

@ -3,6 +3,7 @@
set(ALL_PACKAGES
ADIOS
AMOEBA
ASPHERE
ATC
AWPMD

1
cmake/presets/all_on.cmake
View File

@ -5,6 +5,7 @@
set(ALL_PACKAGES
ADIOS
AMOEBA
ASPHERE
ATC
AWPMD

1
cmake/presets/mingw-cross.cmake
View File

@ -1,4 +1,5 @@
set(WIN_PACKAGES
AMOEBA
ASPHERE
ATC
AWPMD

1
cmake/presets/most.cmake
View File

@ -3,6 +3,7 @@
# are removed. The resulting binary should be able to run most inputs.
set(ALL_PACKAGES
AMOEBA
ASPHERE
BOCS
BODY

1
cmake/presets/windows.cmake
View File

@ -1,4 +1,5 @@
set(WIN_PACKAGES
AMOEBA
ASPHERE
BOCS
BODY

22
doc/src/Howto_amoeba.rst
View File

@ -34,9 +34,9 @@ details for HIPPO are in this paper: :ref:`(Rackers)
.. math::
U & = U_{intermolecular} + U_{intramolecular) \\
U_{intermolecular} & = U_{hal} + U_{repulsion} + U_{dispersion} + U_{multipole} + U_{polar} + U_{qxfer} \\
U_{intramolecular} & = U_{bond} + U_{angle} + U_{torsion} + U_{oop} + U_{b\theta} + U_{UB} + U_{pitorsion} + U_{bitorsion}
U & = U_{intermolecular} + U_{intramolecular} \\
U_{intermolecular} & = U_{hal} + U_{repulsion} + U_{dispersion} + U_{multipole} + U_{polar} + U_{qxfer} \\
U_{intramolecular} & = U_{bond} + U_{angle} + U_{torsion} + U_{oop} + U_{b\theta} + U_{UB} + U_{pitorsion} + U_{bitorsion}
For intermolecular terms, the AMOEBA force field includes only the
:math:`U_{hal}`, :math:`U_{multipole}`, :math:`U_{polar}` terms. The
@ -55,12 +55,12 @@ for a Urey-Bradley bond contribution between the I,K atoms in the IJK
angle.
The :math:`U_{pitorsion}` term is computed by the :doc:`fix
amoeba/pitorsion <fix_pitorsion>` command. It computes 6-body
amoeba/pitorsion <fix_amoeba_pitorsion>` command. It computes 6-body
interaction between a pair of bonded atoms which each have 2
additional bond partners.
The :math:`U_{bitorsion}` term is computed by the :doc:`fix
amoeba/bitorsion <fix_bitorsion>` command. It computes 5-body
amoeba/bitorsion <fix_amoeba_bitorsion>` command. It computes 5-body
interaction between two 4-body torsions (dihedrals) which overlap,
having 3 atoms in common.
@ -72,8 +72,8 @@ compute:
* :doc:`angle_style amoeba <angle_amoeba>`
* :doc:`dihedral_style fourier <dihedral_fourier>`
* :doc:`improper_style amoeba <improper_amoeba>`
* :doc:`fix amoeba/pitorsion <fix_pitorsion>`
* :doc:`fix amoeba/bitorsion <fix_bitorsion>`
* :doc:`fix amoeba/pitorsion <fix_amoeba_pitorsion>`
* :doc:`fix amoeba/bitorsion <fix_amoeba_bitorsion>`
----------
@ -94,14 +94,14 @@ and HIPPO.
improper_style amoeba
# required per-atom data
fix amtype all property/atom i_amtype ghost yes
fix extra all property/atom &
fix extra all property/atom &
i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes
fix polaxe all property/atom i_polaxe
fix pit all amoeba/pitorsion # PiTorsion terms in FF
fix_modify pit energy yes
# Bitorsion terms in FF
fix bit all amoeba/bitorsion bitorsion.ubiquitin.data
fix bit all amoeba/bitorsion bitorsion.ubiquitin.data
fix_modify bit energy yes
read_data data.ubiquitin fix amtype NULL "Tinker Types" &
@ -127,7 +127,7 @@ In the example above the fix ID is *amtype*.
Similarly, if the system you are simulating defines AMOEBA/HIPPO
pitorsion or bitorsion interactions, there will be entries in the data
file for those interactions. They require a :doc:`fix
amoeba/pitortion <fix_amoeba_pitortion>` and :doc:`fix
amoeba/pitortion <fix_amoeba_pitorsion>` and :doc:`fix
amoeba/bitorsion <fix_amoeba_bitorsion>` command be defined. In the
example above, the IDs for these two fixes are *pit* and *bit*.
@ -263,7 +263,7 @@ Switches and their arguments may be specified in any order.
The -xyz switch is required and specifies an input XYZ file as an
argument. The format of this file is an extended XYZ format used by
Tinker for its input. Example *.xyz files are in the examples/amoeba
Tinker for its input. Example \*.xyz files are in the examples/amoeba
directory. The file lists the atoms in the system. Each atom has the
following information: Tinker species name (ignored by LAMMPS), xyz
coordinates, Tinker numeric type, and a list of atom IDs the atom is

19
doc/src/Packages_details.rst
View File

@ -27,6 +27,7 @@ page gives those details.
:columns: 6
* :ref:`ADIOS <PKG-ADIOS>`
* :ref:`AMOEBA <PKG-AMOEBA>`
* :ref:`ASPHERE <PKG-ASPHERE>`
* :ref:`ATC <PKG-ATC>`
* :ref:`AWPMD <PKG-AWPMD>`
@ -148,6 +149,24 @@ This package has :ref:`specific installation instructions <adios>` on the :doc:`
----------
.. _PKG-AMOEBA:
AMOEBA package
---------------
**Contents:**
TODO
**Supporting info:**
* src/AMOEBA: filenames -> commands
* :doc:`AMOEBA and HIPPO howto <Howto_amoeba>`
* examples/amoeba
* TODO
----------
.. _PKG-ASPHERE:
ASPHERE package

5
doc/src/Packages_list.rst
View File

@ -33,6 +33,11 @@ whether an extra library is needed to build and use the package:
- :doc:`dump adios <dump_adios>`
- PACKAGES/adios
- ext
* - :ref:`AMOEBA <PKG-AMOEBA>`
- AMOEBA and HIPPO force fields
- :doc:`AMOEBA and HIPPO howto <Howto_amoeba>`
- amoeba
- no
* - :ref:`ASPHERE <PKG-ASPHERE>`
- aspherical particle models
- :doc:`Howto spherical <Howto_spherical>`

8
doc/src/angle_amoeba.rst
View File

@ -29,9 +29,9 @@ The *amoeba* angle style uses the potential
.. math::
E & = E_a + E_{ba} + E_{ub} \\
E_a = K_2\left(\theta - \theta_0\right)^2 + K_3\left(\theta - \theta_0\right)^3 + K_4\left(\theta - \theta_0\right)^4 + K_5\left(\theta - \theta_0\right)^5 + K_6\left(\theta - \theta_0\right)^6
E_{ba} & = N_1 (r_{ij} - r_1) (\theta - \theta_0) + N_2(r_{jk} - r_2)(\theta - \theta_0)
E_{UB} & = K_{ub} (r_{ik} - r_{ub})^2)
E_a & = K_2\left(\theta - \theta_0\right)^2 + K_3\left(\theta - \theta_0\right)^3 + K_4\left(\theta - \theta_0\right)^4 + K_5\left(\theta - \theta_0\right)^5 + K_6\left(\theta - \theta_0\right)^6 \\
E_{ba} & = N_1 (r_{ij} - r_1) (\theta - \theta_0) + N_2(r_{jk} - r_2)(\theta - \theta_0) \\
E_{UB} & = K_{ub} (r_{ik} - r_{ub})^2
where :math:`E_a` is the angle term, :math:`E_{ba}` is a bond-angle
term, :math:`E_{UB}` is a Urey-Bradley bond term, :math:`\theta_0` is
@ -45,7 +45,7 @@ calculations for the AMOEBA and HIPPO force fields. See the
the implemention of AMOEBA and HIPPO in LAMMPS.
Note that the :math:`E_a` and :math:`E_{ba}` formulas are identical to
those used for the :doc:`angle_style class2/p6 <angle_class2_p6>`
those used for the :doc:`angle_style class2/p6 <angle_class2>`
command, however there is no bond-bond cross term formula for
:math:`E_{bb}`. Additionally, there is a :math:`E_{UB}` term for a
Urey-Bradley bond. It is effectively a harmonic bond between the I

2
doc/src/fix_amoeba_bitorsion.rst
View File

@ -1,7 +1,7 @@
.. index:: fix amoeba/bitorsion
fix amoeba/bitorsion command
================
============================
Syntax
""""""

4
doc/src/fix_amoeba_pitorsion.rst
View File

@ -1,7 +1,7 @@
.. index:: fix amoeba/pitorsion
fix amoeba/pitorsion command
================
============================
Syntax
""""""
@ -86,7 +86,7 @@ the pitorsion 5-atom tuples; it is ignored by LAMMPS. The second
column is the "type" of the interaction; it is an index into the
PiTorsion Coeffs. The remaining 5 columns are the atom IDs of the
atoms in the two 4-atom dihedrals that overlap to create the pitorsion
5-body interaction.
5-body interaction.
Note that the *pitorsion types* and *pitorsions* and *PiTorsion
Coeffs* and *PiTorsions* keywords for the header and body sections of

33
doc/src/pair_amoeba.rst
View File

@ -5,7 +5,7 @@ pair_style amoeba command
=========================
pair_style hippo command
=========================
========================
Syntax
""""""
@ -36,8 +36,8 @@ Additional info
* :doc:`Howto amoeba <howto_ameoba>`
* examples/amoeba
* tools/amoeba
* potentials/*.amoeba
* potentials/*.hippo
* potentials/\*.amoeba
* potentials/\*.hippo
Description
"""""""""""
@ -72,7 +72,7 @@ while the HIPPO force field contains these terms:
U_{hippo} = U_{repulsion} + U_{dispersion} + U_{multipole} + U_{polar} + U_{qxfer}
Conceptually, these terms compute the following interactions:
* :math:`U_{hal}` = buffered 14-7 van der Waals with offsets applied to hydrogen atoms
* :math:`U_{repulsion}` = Pauli repulsion due to rearrangement of electron density
* :math:`U_{dispersion}` = dispersion between correlated, instantaneous induced dipole moments
@ -114,12 +114,12 @@ Only a single pair_coeff command is used with either the *amoeba* and
pair_coeff * * ../potentials/water.prm.hippo ../potentials/water.key.hippo
Examples of the PRM files are in the potentials directory with an
*.amoeba or *.hippo suffix. The examples/amoeba directory has
\*.amoeba or \*.hippo suffix. The examples/amoeba directory has
examples of both PRM and KEY files.
A Tinker PRM file is composed of sections, each of which has multiple
lines. A Tinker KEY file is composed of lines, each of which has a
keyword, followed by zero or more parameters.
keyword followed by zero or more parameters.
The list of PRM sections and KEY keywords which LAMMPS recognizes are
listed on the :doc:`Howto amoeba <howto_ameoba>` doc page. If not
@ -165,8 +165,8 @@ can only use these pair styles with real units.
These potentials do not yet calculate per-atom energy or virial
contributions.
As explained on the :doc:`Howto amoeba <howto_ameoba>` doc page, use
of these pair styles to run a simulation with the AMOEBA or HIPPO
As explained on the :doc:`AMOEBA and HIPPO howto <Howto_amoeba>` page,
use of these pair styles to run a simulation with the AMOEBA or HIPPO
force fields requires several things.
The first is a data file generated by the tools/tinker/tinker2lmp.py
@ -186,18 +186,23 @@ commands for intramolecular interactions may also be required:
* :doc:`angle_style amoeba <angle_amoeba>`
* :doc:`dihedral_style fourier <dihedral_fourier>`
* :doc:`improper_style amoeba <improper_amoeba>`
* :doc:`fix amoeba/pitorsion <fix_pitorsion>`
* :doc:`fix amoeba/bitorsion <fix_bitorsion>`
* :doc:`fix amoeba/pitorsion <fix_amoeba_pitorsion>`
* :doc:`fix amoeba/bitorsion <fix_amoeba_bitorsion>`
----------
Related commands
""""""""""""""""
:doc:`atom_style amoeba <atom_style>`, `bond_style amoeba
:doc:<bond_amoeba>`, `angle_style amoeba <angle_amoeba>`,
:doc:`dihedral_style amoeba <dihedral_amoeba>`, `special_bonds
:doc:one/five <special_bonds>`
:doc:`atom_style amoeba <atom_style>`,
:doc:`bond_style class2 <bond_class2>`,
:doc:`angle_style amoeba <angle_amoeba>`,
:doc:`dihedral_style fourier <dihedral_fourier>`,
:doc:`improper_style amoeba <improper_amoeba>`,
:doc:`fix amoeba/pitorsion <fix_amoeba_pitorsion>`,
:doc:`fix amoeba/bitorsion <fix_amoeba_bitorsion>`,
:doc:`special_bonds one/five <special_bonds>`,
:doc:`fix property/atom <fix_property_atom>`
Default
"""""""

2
doc/src/pair_style.rst
View File

@ -116,6 +116,7 @@ accelerated styles exist.
* :doc:`agni <pair_agni>` - AGNI machine-learning potential
* :doc:`airebo <pair_airebo>` - AIREBO potential of Stuart
* :doc:`airebo/morse <pair_airebo>` - AIREBO with Morse instead of LJ
* :doc:`amoeba <pair_amoeba>` -
* :doc:`atm <pair_atm>` - Axilrod-Teller-Muto potential
* :doc:`awpmd/cut <pair_awpmd>` - Antisymmetrized Wave Packet MD potential for atoms and electrons
* :doc:`beck <pair_beck>` - Beck potential
@ -202,6 +203,7 @@ accelerated styles exist.
* :doc:`hbond/dreiding/lj <pair_hbond_dreiding>` - DREIDING hydrogen bonding LJ potential
* :doc:`hbond/dreiding/morse <pair_hbond_dreiding>` - DREIDING hydrogen bonding Morse potential
* :doc:`hdnnp <pair_hdnnp>` - High-dimensional neural network potential
* :doc:`hippo <pair_amoeba>` -
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` - registry-dependent interlayer potential (ILP)
* :doc:`ilp/tmd <pair_ilp_tmd>` - interlayer potential (ILP) potential for transition metal dichalcogenides (TMD)
* :doc:`kim <pair_kim>` - interface to potentials provided by KIM project

6
src/AMOEBA/amoeba_charge_transfer.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -118,7 +118,7 @@ void PairAmoeba::charge_transfer()
de = felec * de * factor_mpole;
// use energy switching if near the cutoff distance
if (r2 > cut2) {
r3 = r2 * r;
r4 = r2 * r2;
@ -165,7 +165,7 @@ void PairAmoeba::charge_transfer()
// energy = e
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}

262
src/AMOEBA/amoeba_convolution.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -58,8 +58,8 @@ enum{MPOLE_GRID,POLAR_GRID,POLAR_GRIDC,DISP_GRID,INDUCE_GRID,INDUCE_GRIDC};
------------------------------------------------------------------------- */
AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
int nx_caller, int ny_caller, int nz_caller,
int order_caller, int which_caller) :
int nx_caller, int ny_caller, int nz_caller,
int order_caller, int which_caller) :
Pointers(lmp)
{
amoeba = pair;
@ -73,9 +73,9 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
if (which == POLAR_GRIDC || which == INDUCE_GRIDC) flag3d = 0;
// NOTE: worry about overflow
nfft_global = nx * ny * nz;
// global indices of grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global grid that I own without ghost cells
@ -112,10 +112,10 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
// convert to triclinic if necessary
// nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping
// NOTE: this needs to be computed same as IGRID in amoeba
double *prd,*boxlo,*sublo,*subhi;
int triclinic = domain->triclinic;
if (triclinic == 0) {
prd = domain->prd;
boxlo = domain->boxlo;
@ -153,7 +153,7 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
nhi = static_cast<int> ((subhi[2]+dist[2]-boxlo[2]) * nz/zprd);
nzlo_out = nlo + nlower;
nzhi_out = nhi + nupper;
// x-pencil decomposition of FFT mesh
// global indices range from 0 to N-1
// each proc owns entire x-dimension, clumps of columns in y,z dimensions
@ -166,7 +166,7 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
int me = comm->me;
int nprocs = comm->nprocs;
int npey_fft,npez_fft;
if (nz >= nprocs) {
npey_fft = 1;
@ -189,7 +189,7 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
// nfft_owned = owned grid points in FFT decomp
// ngrid_either = max of nbrick_onwed and nfft_owned
// nfft = total FFT grid points
nbrick_owned = (nxhi_in-nxlo_in+1) * (nyhi_in-nylo_in+1) *
(nzhi_in-nzlo_in+1);
nbrick_ghosts = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
@ -202,43 +202,43 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair,
// instantiate FFT, GridComm, and Remap
int tmp;
fft1 = new FFT3d(lmp,world,nx,ny,nz,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
1,0,&tmp,0);
// 0,0,&tmp,0);
// 0,0,&tmp,0);
fft2 = new FFT3d(lmp,world,nx,ny,nz,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
//1,0,&tmp,0);
0,0,&tmp,0);
gc = new GridComm(lmp,world,nx,ny,nz,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int nqty = flag3d ? 1 : 2;
remap = new Remap(lmp,world,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nqty,0,0,FFT_PRECISION,0);
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
nqty,0,0,FFT_PRECISION,0);
// memory allocations
if (flag3d) {
memory->create3d_offset(grid_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"amoeba:grid_brick");
nxlo_out,nxhi_out,"amoeba:grid_brick");
grid_brick_start = &grid_brick[nzlo_out][nylo_out][nxlo_out];
cgrid_brick = NULL;
} else {
memory->create4d_offset_last(cgrid_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,2,"amoeba:cgrid_brick");
nxlo_out,nxhi_out,2,"amoeba:cgrid_brick");
grid_brick_start = &cgrid_brick[nzlo_out][nylo_out][nxlo_out][0];
grid_brick = NULL;
}
memory->create(grid_fft,ngrid_either,"amoeba:grid_fft");
memory->create(cfft,2*ngrid_either,"amoeba:cfft");
@ -263,7 +263,7 @@ AmoebaConvolution::~AmoebaConvolution()
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
memory->destroy(remap_buf);
delete fft1;
delete fft2;
delete gc;
@ -288,7 +288,7 @@ void *AmoebaConvolution::zero_3d()
{
if (!grid_brick) return NULL;
memset(&(grid_brick[nzlo_out][nylo_out][nxlo_out]),0,
nbrick_ghosts*sizeof(FFT_SCALAR));
nbrick_ghosts*sizeof(FFT_SCALAR));
return (void *) grid_brick;
}
@ -298,7 +298,7 @@ void *AmoebaConvolution::zero_4d()
{
if (!cgrid_brick) return NULL;
memset(&(cgrid_brick[nzlo_out][nylo_out][nxlo_out][0]),0,
2*nbrick_ghosts*sizeof(FFT_SCALAR));
2*nbrick_ghosts*sizeof(FFT_SCALAR));
return (void *) cgrid_brick;
}
@ -327,11 +327,11 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_3d()
if (DEBUG) debug_scalar(GRIDBRICK_OUT,"PRE Convo / PRE GridComm");
gc->reverse_comm(GridComm::PAIR,amoeba,1,sizeof(FFT_SCALAR),which,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
if (DEBUG) debug_scalar(GRIDBRICK_IN,"PRE Convo / POST GridComm");
if (DEBUG) debug_file(GRIDBRICK_IN,"pre.convo.post.gridcomm");
// copy owned 3d brick grid values to FFT grid
n = 0;
@ -339,7 +339,7 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_3d()
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
grid_fft[n++] = grid_brick[iz][iy][ix];
// remap FFT grid from brick to x pencil partitioning
remap->perform(grid_fft,grid_fft,remap_buf);
@ -348,7 +348,7 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_3d()
if (DEBUG) debug_file(FFT,"pre.convo.post.remap");
// copy real values into complex grid
n = 0;
for (int i = 0; i < nfft_owned; i++) {
cfft[n++] = grid_fft[i];
@ -356,7 +356,7 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_3d()
}
// perform forward FFT
fft1->compute(cfft,cfft,FFT3d::FORWARD);
if (SCALE) {
@ -383,7 +383,7 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_4d()
if (DEBUG) debug_scalar(GRIDBRICK_OUT,"PRE Convo / PRE GridComm");
gc->reverse_comm(GridComm::PAIR,amoeba,2,sizeof(FFT_SCALAR),which,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
if (DEBUG) debug_scalar(GRIDBRICK_IN,"PRE Convo / POST GridComm");
if (DEBUG) debug_file(GRIDBRICK_IN,"pre.convo.post.gridcomm");
@ -400,14 +400,14 @@ FFT_SCALAR *AmoebaConvolution::pre_convolution_4d()
// remap FFT grid from brick to x pencil partitioning
// NOTE: could just setup FFT to start from brick decomp and skip remap
remap->perform(cfft,cfft,remap_buf);
if (DEBUG) debug_scalar(FFT,"PRE Convo / POST Remap");
if (DEBUG) debug_file(FFT,"pre.convo.post.remap");
// perform forward FFT
fft1->compute(cfft,cfft,FFT3d::FORWARD);
if (SCALE) {
@ -440,7 +440,7 @@ void *AmoebaConvolution::post_convolution()
void *AmoebaConvolution::post_convolution_3d()
{
int ix,iy,iz,n;
// perform backward FFT
if (DEBUG) debug_scalar(CFFT1,"POST Convo / PRE FFT");
@ -457,17 +457,17 @@ void *AmoebaConvolution::post_convolution_3d()
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
grid_brick[iz][iy][ix] = cfft[n];
n += 2;
grid_brick[iz][iy][ix] = cfft[n];
n += 2;
}
// forward comm to populate ghost grid values
if (DEBUG) debug_scalar(GRIDBRICK_IN,"POST Convo / PRE gridcomm");
if (DEBUG) debug_file(GRIDBRICK_IN,"post.convo.pre.gridcomm");
gc->forward_comm(GridComm::PAIR,amoeba,1,sizeof(FFT_SCALAR),which,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
return (void *) grid_brick;
}
@ -496,17 +496,17 @@ void *AmoebaConvolution::post_convolution_4d()
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
cgrid_brick[iz][iy][ix][0] = cfft[n++];
cgrid_brick[iz][iy][ix][1] = cfft[n++];
cgrid_brick[iz][iy][ix][0] = cfft[n++];
cgrid_brick[iz][iy][ix][1] = cfft[n++];
}
// forward comm to populate ghost grid values
if (DEBUG) debug_scalar(GRIDBRICK_IN,"POST Convo / PRE gridcomm");
if (DEBUG) debug_file(GRIDBRICK_IN,"post.convo.pre.gridcomm");
gc->forward_comm(GridComm::PAIR,amoeba,2,sizeof(FFT_SCALAR),which,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
return (void *) cgrid_brick;
}
@ -524,7 +524,7 @@ void AmoebaConvolution::kspacebbox(double r, double *b)
{
double *h = domain->h;
double lx,ly,lz,xy,xz,yz;
lx = h[0]; ly = h[1]; lz = h[2];
yz = h[3]; xz = h[4]; xy = h[5];
@ -585,42 +585,42 @@ void AmoebaConvolution::debug_scalar(int array, const char *label)
if (array == GRIDBRICK_OUT) {
if (flag3d) {
for (int iz = nzlo_out; iz <= nzhi_out; iz++)
for (int iy = nylo_out; iy <= nyhi_out; iy++)
for (int ix = nxlo_out; ix <= nxhi_out; ix++)
sum += grid_brick[iz][iy][ix];
for (int iy = nylo_out; iy <= nyhi_out; iy++)
for (int ix = nxlo_out; ix <= nxhi_out; ix++)
sum += grid_brick[iz][iy][ix];
} else {
for (int iz = nzlo_out; iz <= nzhi_out; iz++)
for (int iy = nylo_out; iy <= nyhi_out; iy++)
for (int ix = nxlo_out; ix <= nxhi_out; ix++) {
sum += cgrid_brick[iz][iy][ix][0];
sum += cgrid_brick[iz][iy][ix][1];
}
for (int iy = nylo_out; iy <= nyhi_out; iy++)
for (int ix = nxlo_out; ix <= nxhi_out; ix++) {
sum += cgrid_brick[iz][iy][ix][0];
sum += cgrid_brick[iz][iy][ix][1];
}
}
}
if (array == GRIDBRICK_IN) {
if (flag3d) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
sum += grid_brick[iz][iy][ix];
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
sum += grid_brick[iz][iy][ix];
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
sum += cgrid_brick[iz][iy][ix][0];
sum += cgrid_brick[iz][iy][ix][1];
}
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
sum += cgrid_brick[iz][iy][ix][0];
sum += cgrid_brick[iz][iy][ix][1];
}
}
}
if (array == FFT) {
if (flag3d) {
for (int i = 0; i < nfft_owned; i++)
sum += grid_fft[i];
sum += grid_fft[i];
} else {
for (int i = 0; i < 2*nfft_owned; i++)
sum += cfft[i];
sum += cfft[i];
}
}
@ -628,13 +628,13 @@ void AmoebaConvolution::debug_scalar(int array, const char *label)
for (int i = 0; i < 2*nfft_owned; i++)
sum += cfft[i];
}
if (array == CFFT2) {
for (int i = 0; i < 2*nbrick_owned; i++)
sum += cfft[i];
}
/*
/*
double sumall;
MPI_Allreduce(&sum,&sumall,1,MPI_DOUBLE,MPI_SUM,world);
if (comm->me == 0) printf("%s: %s: %12.8g\n",labels[which],label,sumall);
@ -654,14 +654,14 @@ void AmoebaConvolution::debug_file(int array, const char *label)
int nprocs = comm->nprocs;
// open file
char fname[128];
sprintf(fname,"tmp.%s.%s",labels[which],label);
if (me == 0) fp = fopen(fname,"w");
// file header
// ncol = # of columns, including grid cell ID
bigint ntot = nx * ny * nz;
int ncol;
@ -701,58 +701,58 @@ void AmoebaConvolution::debug_file(int array, const char *label)
int ngridmax;
MPI_Allreduce(&ngrid,&ngridmax,1,MPI_INT,MPI_MAX,world);
double *buf,*buf2;
memory->create(buf,ncol*ngridmax,"amoeba:buf");
memory->create(buf2,ncol*ngridmax,"amoeba:buf2");
ngrid = 0;
if (array == GRIDBRICK_IN) {
if (flag3d) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = grid_brick[iz][iy][ix];
ngrid++;
}
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = grid_brick[iz][iy][ix];
ngrid++;
}
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cgrid_brick[iz][iy][ix][0];
buf[ncol*ngrid+2] = cgrid_brick[iz][iy][ix][1];
ngrid++;
}
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cgrid_brick[iz][iy][ix][0];
buf[ncol*ngrid+2] = cgrid_brick[iz][iy][ix][1];
ngrid++;
}
}
}
if (array == FFT) {
if (flag3d) {
int m = 0;
for (int iz = nzlo_fft; iz <= nzhi_fft; iz++)
for (int iy = nylo_fft; iy <= nyhi_fft; iy++)
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = grid_fft[m++];
ngrid++;
}
} else {
for (int iy = nylo_fft; iy <= nyhi_fft; iy++)
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = grid_fft[m++];
ngrid++;
}
} else {
int m = 0;
for (int iz = nzlo_fft; iz <= nzhi_fft; iz++)
for (int iy = nylo_fft; iy <= nyhi_fft; iy++)
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
for (int iy = nylo_fft; iy <= nyhi_fft; iy++)
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
}
}
@ -760,62 +760,62 @@ void AmoebaConvolution::debug_file(int array, const char *label)
int m = 0;
for (int iz = nzlo_fft; iz <= nzhi_fft; iz++)
for (int iy = nylo_fft; iy <= nyhi_fft; iy++)
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
for (int ix = nxlo_fft; ix <= nxhi_fft; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
}
if (array == CFFT2) {
int m = 0;
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
int id = iz*ny*nx + iy*nx + ix + 1;
buf[ncol*ngrid] = id;
buf[ncol*ngrid+1] = cfft[m++];
buf[ncol*ngrid+2] = cfft[m++];
ngrid++;
}
}
// proc 0 outputs values
// pings other procs, send/recv of their values
int tmp,nlines;
MPI_Request request;
MPI_Status status;
if (me == 0) {
for (int iproc = 0; iproc < nprocs; iproc++) {
if (iproc) {
MPI_Irecv(buf,ngridmax*ncol,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,me+iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= ncol;
MPI_Irecv(buf,ngridmax*ncol,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,me+iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= ncol;
} else nlines = ngrid;
int n = 0;
for (int m = 0; m < nlines; m++) {
if (ncol == 2)
fprintf(fp,"%d %12.8g\n",(int) buf[n],buf[n+1]);
else if (ncol == 3)
fprintf(fp,"%d %12.8g %12.8g\n",(int ) buf[n],buf[n+1],buf[n+2]);
n += ncol;
if (ncol == 2)
fprintf(fp,"%d %12.8g\n",(int) buf[n],buf[n+1]);
else if (ncol == 3)
fprintf(fp,"%d %12.8g %12.8g\n",(int ) buf[n],buf[n+1],buf[n+2]);
n += ncol;
}
}
} else {
MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(buf,ngrid*ncol,MPI_DOUBLE,0,0,world);
}
// close file
if (me == 0) fclose(fp);
// clean up

14
src/AMOEBA/amoeba_convolution.h
View File

@ -1,6 +1,6 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -39,14 +39,14 @@ class AmoebaConvolution : protected Pointers {
int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
int nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft;
double *grid_brick_start; // lower left corner of (c)grid_brick data
AmoebaConvolution(class LAMMPS *, class Pair *,
int, int, int, int, int);
int, int, int, int, int);
~AmoebaConvolution();
void *zero();
FFT_SCALAR *pre_convolution();
void *post_convolution();
private:
int which; // caller name for convolution being performed
int flag3d; // 1 if using 3d grid_brick, 0 for 4d cgrid_brick
@ -58,13 +58,13 @@ class AmoebaConvolution : protected Pointers {
class FFT3d *fft1,*fft2;
class GridComm *gc;
class Remap *remap;
double ***grid_brick; // 3d real brick grid with ghosts
double ****cgrid_brick; // 4d complex brick grid with ghosts
FFT_SCALAR *grid_fft; // 3d FFT grid as 1d vector
FFT_SCALAR *cfft; // 3d complex FFT grid as 1d vector
double *gc_buf1,*gc_buf2; // buffers for GridComm
double *remap_buf; // buffer for Remap

98
src/AMOEBA/amoeba_dispersion.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -97,7 +97,7 @@ void PairAmoeba::dispersion_real()
int *ilist,*jlist,*numneigh,**firstneigh;
// owned atoms
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
@ -172,10 +172,10 @@ void PairAmoeba::dispersion_real()
tk2 = tk * tk;
damp3 = 1.0 - ti2*(1.0+di+0.5*di2)*expi - tk2*(1.0+dk+0.5*dk2)*expk -
2.0*ti2*tk*(1.0+di)*expi - 2.0*tk2*ti*(1.0+dk)*expk;
damp5 = 1.0 - ti2*(1.0+di+0.5*di2+di3/6.0)*expi -
tk2*(1.0+dk+0.5*dk2 + dk3/6.0)*expk -
damp5 = 1.0 - ti2*(1.0+di+0.5*di2+di3/6.0)*expi -
tk2*(1.0+dk+0.5*dk2 + dk3/6.0)*expk -
2.0*ti2*tk*(1.0+di+di2/3.0)*expi - 2.0*tk2*ti*(1.0+dk+dk2/3.0)*expk;
ddamp = 0.25 * di2 * ti2 * ai * expi * (r*ai+4.0*tk-1.0) +
ddamp = 0.25 * di2 * ti2 * ai * expi * (r*ai+4.0*tk-1.0) +
0.25 * dk2 * tk2 * ak * expk * (r*ak+4.0*ti-1.0);
} else {
@ -187,7 +187,7 @@ void PairAmoeba::dispersion_real()
}
damp = 1.5*damp5 - 0.5*damp3;
// apply damping and scaling factors for this interaction
scale = factor_disp * damp*damp;
@ -229,7 +229,7 @@ void PairAmoeba::dispersion_real()
// energy = e
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}
@ -298,11 +298,11 @@ void PairAmoeba::dispersion_kspace()
// pre-convolution operations including forward FFT
// gridfft = my portion of complex 3d grid in FFT decomposition
double *gridfft = d_kspace->pre_convolution();
// ---------------------
// convolution operation
// convolution operation
// ---------------------
nhalf1 = (nfft1+1) / 2;
@ -329,39 +329,39 @@ void PairAmoeba::dispersion_kspace()
for (k = nzlo; k <= nzhi; k++) {
for (j = nylo; j <= nyhi; j++) {
for (i = nxlo; i <= nxhi; i++) {
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
h = sqrt(hsq);
b = h*bfac;
hhh = h*hsq;
term = -b*b;
expterm = 0.0;
erfcterm = erfc(b);
denom = denom0*bsmod1[i]*bsmod2[j]*bsmod3[k];
if (term > -50.0 && hsq != 0.0) {
expterm = exp(term);
erfcterm = erfc(b);
term1 = fac1*erfcterm*hhh + expterm*(fac2 + fac3*hsq);
struc2 = gridfft[n]*gridfft[n] + gridfft[n+1]*gridfft[n+1];
e = -(term1 / denom) * struc2;
edisp += e;
vterm = 3.0 * (fac1*erfcterm*h + fac3*expterm) * struc2/denom;
virdisp[0] += h1*h1*vterm - e;
virdisp[1] += h2*h2*vterm - e;
virdisp[2] += h3*h3*vterm - e;
virdisp[3] += h1*h2*vterm;
virdisp[4] += h1*h3*vterm;
virdisp[5] += h2*h3*vterm;
} else term1 = 0.0;
// NOTE: pre-calc this division only once
gridfft[n] *= -(term1/denom);
gridfft[n+1] *= -(term1/denom);
n += 2;
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
h = sqrt(hsq);
b = h*bfac;
hhh = h*hsq;
term = -b*b;
expterm = 0.0;
erfcterm = erfc(b);
denom = denom0*bsmod1[i]*bsmod2[j]*bsmod3[k];
if (term > -50.0 && hsq != 0.0) {
expterm = exp(term);
erfcterm = erfc(b);
term1 = fac1*erfcterm*hhh + expterm*(fac2 + fac3*hsq);
struc2 = gridfft[n]*gridfft[n] + gridfft[n+1]*gridfft[n+1];
e = -(term1 / denom) * struc2;
edisp += e;
vterm = 3.0 * (fac1*erfcterm*h + fac3*expterm) * struc2/denom;
virdisp[0] += h1*h1*vterm - e;
virdisp[1] += h2*h2*vterm - e;
virdisp[2] += h3*h3*vterm - e;
virdisp[3] += h1*h2*vterm;
virdisp[4] += h1*h3*vterm;
virdisp[5] += h2*h3*vterm;
} else term1 = 0.0;
// NOTE: pre-calc this division only once
gridfft[n] *= -(term1/denom);
gridfft[n+1] *= -(term1/denom);
n += 2;
}
}
}
@ -371,7 +371,7 @@ void PairAmoeba::dispersion_kspace()
double ***gridpost = (double ***) d_kspace->post_convolution();
// get first derivatives of the reciprocal space energy
// get first derivatives of the reciprocal space energy
int nlpts = (bsorder-1) / 2;
int nrpts = bsorder - nlpts - 1;
@ -391,17 +391,17 @@ void PairAmoeba::dispersion_kspace()
t2 = thetai2[m][jb][0];
dt2 = nfft2 * thetai2[m][jb][1];
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t1 = thetai1[m][ib][0];
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t1 = thetai1[m][ib][0];
dt1 = nfft1 * thetai1[m][ib][1];
term = gridpost[k][j][i];
de1 += 2.0*term*dt1*t2*t3;
de2 += 2.0*term*dt2*t1*t3;
de3 += 2.0*term*dt3*t1*t2;
i++;
i++;
}
j++;
j++;
}
k++;
}
@ -415,7 +415,7 @@ void PairAmoeba::dispersion_kspace()
// account for the energy and virial correction terms
term = csixpr * aewald*aewald*aewald / denom0;
if (me == 0) {
edisp -= term;
virdisp[0] += term;

136
src/AMOEBA/amoeba_file.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel ator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -135,7 +135,7 @@ void PairAmoeba::read_prmfile(char *filename)
if (forcefield_flag == 0 && section != FFIELD)
error->all(FLERR,"Force Field is not first section of "
"pair amoeba potential file");
if (section != FFIELD && section != LITERATURE && section != ATOMTYPE &&
if (section != FFIELD && section != LITERATURE && section != ATOMTYPE &&
section != UNKNOWN && atomtype_flag == 0)
error->all(FLERR,"Atom Type section of pair amoeba potential file "
"must come before all but Force Field section");
@ -159,8 +159,8 @@ void PairAmoeba::read_prmfile(char *filename)
// convert all chars in line to lower-case
for (int i = 0; i < n; i++)
line[i] = tolower(line[i]);
line[i] = tolower(line[i]);
char *copy;
char **words;
int nwords = tokenize(line,words,copy);
@ -205,11 +205,11 @@ void PairAmoeba::read_prmfile(char *filename)
void PairAmoeba::read_keyfile(char *filename)
{
double aprd,bprd,cprd;
// default settings for which there are keyword options
aprd = bprd = cprd = 0.0;
vdwcut = 9.0;
vdwtaper = 0.9 * vdwcut;
repcut = 6.0;
@ -254,7 +254,7 @@ void PairAmoeba::read_keyfile(char *filename)
// done if keyfile not specified by pair_coeff command
if (!filename) return;
// open key file
int me = comm->me;
@ -320,14 +320,14 @@ void PairAmoeba::read_keyfile(char *filename)
if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid");
double taper = utils::numeric(FLERR,words[1],true,lmp);
if (taper >= 1.0) {
vdwtaper = reptaper = disptaper = mpoletaper = ctrntaper = taper;
vdwtaper = reptaper = disptaper = mpoletaper = ctrntaper = taper;
} else {
taper = -taper;
vdwtaper = taper * vdwcut;
reptaper = taper * repcut;
disptaper = taper * dispcut;
mpoletaper = taper * mpolecut;
ctrntaper = taper * ctrncut;
taper = -taper;
vdwtaper = taper * vdwcut;
reptaper = taper * repcut;
disptaper = taper * dispcut;
mpoletaper = taper * mpolecut;
ctrntaper = taper * ctrncut;
}
} else if (strcmp(keyword,"vdw-cutoff") == 0) {
if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid");
@ -396,7 +396,7 @@ void PairAmoeba::read_keyfile(char *filename)
} else if (strcmp(keyword,"dewald") == 0) {
if (nwords != 1) error->all(FLERR,"AMOEBA keyfile line is invalid");
use_dewald = 1;
} else if (strcmp(keyword,"pme-order") == 0) {
if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid");
bseorder = utils::numeric(FLERR,words[1],true,lmp);
@ -409,24 +409,24 @@ void PairAmoeba::read_keyfile(char *filename)
} else if (strcmp(keyword,"pme-grid") == 0) {
if (nwords != 2 && nwords != 4)
error->all(FLERR,"AMOEBA keyfile line is invalid");
error->all(FLERR,"AMOEBA keyfile line is invalid");
if (nwords == 2)
nefft1 = nefft2 = nefft3 = utils::numeric(FLERR,words[1],true,lmp);
nefft1 = nefft2 = nefft3 = utils::numeric(FLERR,words[1],true,lmp);
else {
nefft1 = utils::numeric(FLERR,words[1],true,lmp);
nefft2 = utils::numeric(FLERR,words[2],true,lmp);
nefft3 = utils::numeric(FLERR,words[3],true,lmp);
nefft1 = utils::numeric(FLERR,words[1],true,lmp);
nefft2 = utils::numeric(FLERR,words[2],true,lmp);
nefft3 = utils::numeric(FLERR,words[3],true,lmp);
}
pmegrid_key = 1;
} else if (strcmp(keyword,"dpme-grid") == 0) {
if (nwords != 2 && nwords != 4)
error->all(FLERR,"AMOEBA keyfile line is invalid");
error->all(FLERR,"AMOEBA keyfile line is invalid");
if (nwords == 2)
ndfft1 = ndfft2 = ndfft3 = utils::numeric(FLERR,words[1],true,lmp);
ndfft1 = ndfft2 = ndfft3 = utils::numeric(FLERR,words[1],true,lmp);
else {
ndfft1 = utils::numeric(FLERR,words[1],true,lmp);
ndfft2 = utils::numeric(FLERR,words[2],true,lmp);
ndfft3 = utils::numeric(FLERR,words[3],true,lmp);
ndfft1 = utils::numeric(FLERR,words[1],true,lmp);
ndfft2 = utils::numeric(FLERR,words[2],true,lmp);
ndfft3 = utils::numeric(FLERR,words[3],true,lmp);
}
dpmegrid_key = 1;
@ -448,27 +448,27 @@ void PairAmoeba::read_keyfile(char *filename)
} else if (strcmp(words[0],"polarization") == 0) {
if (strcmp(words[1],"mutual") == 0) poltyp = MUTUAL;
else if (strstr(words[1],"opt") == words[1]) {
poltyp = OPT;
if (strcmp(words[1],"opt") == 0) optorder = 4;
else optorder = utils::inumeric(FLERR,&words[1][3],true,lmp);
if (optorder < 1 || optorder > 6)
error->all(FLERR,"Unrecognized polarization OPTn in AMOEBA FF file");
poltyp = OPT;
if (strcmp(words[1],"opt") == 0) optorder = 4;
else optorder = utils::inumeric(FLERR,&words[1][3],true,lmp);
if (optorder < 1 || optorder > 6)
error->all(FLERR,"Unrecognized polarization OPTn in AMOEBA FF file");
} else if (strcmp(words[1],"tcg") == 0)
error->all(FLERR,"Polarization TCG not yet supported in AMOEBA/HIPPO");
error->all(FLERR,"Polarization TCG not yet supported in AMOEBA/HIPPO");
else if (strcmp(words[1],"direct") == 0) poltyp = DIRECT;
else error->all(FLERR,"Unrecognized polarization in AMOEBA FF file");
} else if (strcmp(keyword,"polar-predict") == 0) {
if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid");
if (strcmp(words[1],"gear") == 0) {
polpred = GEAR;
maxualt = 7;
polpred = GEAR;
maxualt = 7;
} else if (strcmp(words[1],"aspc") == 0) {
polpred = ASPC;
maxualt = 17;
polpred = ASPC;
maxualt = 17;
} else if (strcmp(words[1],"lsqr") == 0) {
polpred = LSQR;
maxualt = 7;
polpred = LSQR;
maxualt = 7;
} else error->all(FLERR,"AMOEBA keyfile line is invalid");
use_pred = 1;
} else if (strcmp(keyword,"polar-iter") == 0) {
@ -493,7 +493,7 @@ void PairAmoeba::read_keyfile(char *filename)
} else if (strcmp(keyword,"pcg-peek") == 0) {
if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid");
pcgpeek = utils::numeric(FLERR,words[1],true,lmp);
} else {}
delete [] copy;
@ -501,14 +501,14 @@ void PairAmoeba::read_keyfile(char *filename)
}
// close key file
if (me == 0) fclose(fptr);
// cutoff resets for long-range interactions
if (use_ewald) mpolecut = ewaldcut;
if (use_dewald) dispcut = dewaldcut;
// error checks
if (use_ewald || use_dewald) {
@ -555,8 +555,8 @@ int PairAmoeba::read_section_name(FILE *fp, char *line)
}
break;
}
if ((strstr(words[0],"##") != words[0]) ||
if ((strstr(words[0],"##") != words[0]) ||
(strstr(words[nwords-1],"##") != words[nwords-1]))
error->one(FLERR,"Section header of pair amoeba potential file is invalid");
@ -580,7 +580,7 @@ int PairAmoeba::read_section_name(FILE *fp, char *line)
/* ---------------------------------------------------------------------- */
int PairAmoeba::read_section_line(FILE *fp, char *line,
int PairAmoeba::read_section_line(FILE *fp, char *line,
int &nextflag, char *next)
{
// loop on read line
@ -594,7 +594,7 @@ int PairAmoeba::read_section_line(FILE *fp, char *line,
// until a line is blank, starts with #, starts with alpha char, or EOF
// save next line read, and set nextflag for next call
// return length of line
char *ptr,*copy,*copy_next;
char **words,**words_next;
int nwords,nwords_next;
@ -776,7 +776,7 @@ void PairAmoeba::file_ffield(int nwords, char **words)
poltyp = OPT;
if (strcmp(words[1],"opt") == 0) optorder = 4;
else optorder = utils::inumeric(FLERR,&words[1][3],true,lmp);
if (optorder < 1 || optorder > 6)
if (optorder < 1 || optorder > 6)
error->all(FLERR,"Unrecognized polarization OPTn in AMOEBA FF file");
} else if (strcmp(words[1],"tcg") == 0)
error->all(FLERR,"Polarization TCG not yet supported in AMOEBA/HIPPO");
@ -785,7 +785,7 @@ void PairAmoeba::file_ffield(int nwords, char **words)
// NOTE: enable all variants of special settings
// do these need to be set to defaults if don't appear in file?
} else if (strcmp(words[0],"vdw-12-scale") == 0) {
special_hal[1] = utils::numeric(FLERR,words[1],true,lmp);
} else if (strcmp(words[0],"vdw-13-scale") == 0) {
@ -883,9 +883,9 @@ void PairAmoeba::file_atomtype(int nwords, char **words)
{
if (nwords < 8)
error->all(FLERR,"AMOEBA atom type line is invalid");
if (strcmp(words[0],"atom") != 0)
if (strcmp(words[0],"atom") != 0)
error->all(FLERR,"AMOEBA atom type line is invalid");
int itype = utils::inumeric(FLERR,words[1],true,lmp);
int iclass = utils::inumeric(FLERR,words[2],true,lmp);
@ -913,7 +913,7 @@ void PairAmoeba::file_vdwl(int nwords, char **words)
{
if (nwords != 4 && nwords != 5)
error->all(FLERR,"AMOEBA Van der Waals line is invalid");
if (strcmp(words[0],"vdw") != 0)
if (strcmp(words[0],"vdw") != 0)
error->all(FLERR,"AMOEBA Van der Waals line is invalid");
int iclass = utils::inumeric(FLERR,words[1],true,lmp);
@ -930,9 +930,9 @@ void PairAmoeba::file_vdwl(int nwords, char **words)
void PairAmoeba::file_vdwl_pair(int nwords, char **words)
{
if (nwords != 5)
if (nwords != 5)
error->all(FLERR,"AMOEBA Van der Waals pair line is invalid");
if (strcmp(words[0],"vdwpr") != 0)
if (strcmp(words[0],"vdwpr") != 0)
error->all(FLERR,"AMOEBA Van der Waals pair line is invalid");
if (nvdwl_pair == max_vdwl_pair) {
@ -955,7 +955,7 @@ void PairAmoeba::file_bstretch(int nwords, char **words)
{
if (nwords < 5)
error->all(FLERR,"AMOEBA bond stretch line is invalid");
if (strcmp(words[0],"bond") != 0)
if (strcmp(words[0],"bond") != 0)
error->all(FLERR,"AMOEBA bond stretch line is invalid");
}
@ -963,9 +963,9 @@ void PairAmoeba::file_bstretch(int nwords, char **words)
void PairAmoeba::file_sbend(int nwords, char **words)
{
if (nwords != 6)
if (nwords != 6)
error->all(FLERR,"AMOEBA strectch-bend line is invalid");
if (strstr(words[0],"strbnd") != words[0])
if (strstr(words[0],"strbnd") != words[0])
error->all(FLERR,"AMOEBA strectch-bend line is invalid");
}
@ -983,7 +983,7 @@ void PairAmoeba::file_pauli(int nwords, char **words)
{
if (nwords < 5)
error->all(FLERR,"AMOEBA Pauli repulsion line is invalid");
if (strstr(words[0],"repulsion") != words[0])
if (strstr(words[0],"repulsion") != words[0])
error->all(FLERR,"AMOEBA Pauli repulsion line is invalid");
int itype = utils::inumeric(FLERR,words[1],true,lmp);
@ -1003,7 +1003,7 @@ void PairAmoeba::file_dispersion(int nwords, char **words)
{
if (nwords < 4)
error->all(FLERR,"AMOEBA dispersion line is invalid");
if (strstr(words[0],"dispersion") != words[0])
if (strstr(words[0],"dispersion") != words[0])
error->all(FLERR,"AMOEBA dipersion line is invalid");
int iclass = utils::inumeric(FLERR,words[1],true,lmp);
@ -1020,7 +1020,7 @@ void PairAmoeba::file_ub(int nwords, char **words)
{
if (nwords != 6)
error->all(FLERR,"AMOEBA Urey-Bradley line is invalid");
if (strstr(words[0],"ureybrad") != words[0])
if (strstr(words[0],"ureybrad") != words[0])
error->all(FLERR,"AMOEBA Urey-Bradley line is invalid");
}
@ -1030,7 +1030,7 @@ void PairAmoeba::file_outplane(int nwords, char **words)
{
if (nwords != 6)
error->all(FLERR,"AMOEBA out-of-plane bend line is invalid");
if (strstr(words[0],"opbend") != words[0])
if (strstr(words[0],"opbend") != words[0])
error->all(FLERR,"AMOEBA out-of-plane bend line is invalid");
}
@ -1040,7 +1040,7 @@ void PairAmoeba::file_torsion(int nwords, char **words)
{
if (nwords != 14)
error->all(FLERR,"AMOEBA torsional line is invalid");
if (strstr(words[0],"torsion") != words[0])
if (strstr(words[0],"torsion") != words[0])
error->all(FLERR,"AMOEBA torsional line is invalid");
}
@ -1060,7 +1060,7 @@ void PairAmoeba::file_multipole(int nwords, char **words)
{
if (nwords < 12 || nwords > 15)
error->all(FLERR,"AMOEBA atomic multipole line is invalid");
if (strstr(words[0],"multipole") != words[0])
if (strstr(words[0],"multipole") != words[0])
error->all(FLERR,"AMOEBA atomic multipole line is invalid");
int itype = utils::inumeric(FLERR,words[1],true,lmp);
@ -1070,7 +1070,7 @@ void PairAmoeba::file_multipole(int nwords, char **words)
int iframe = nmultiframe[itype];
if (iframe == MAX_FRAME_PER_TYPE)
error->all(FLERR,"AMOEBA MAX_FRAME_PER_TYPE is too small");
int extra;
if (nwords == 12) {
zpole[itype][iframe] = xpole[itype][iframe] = ypole[itype][iframe] = 0;
@ -1116,7 +1116,7 @@ void PairAmoeba::file_multipole(int nwords, char **words)
// quadrupole terms divided by 3 for use as traceless values
for (int i = 1; i < 4; i++)
fpole[itype][iframe][i] *= BOHR;
fpole[itype][iframe][i] *= BOHR;
for (int i = 4; i < 13; i++)
fpole[itype][iframe][i] *= BOHR*BOHR / 3.0;
@ -1135,7 +1135,7 @@ void PairAmoeba::file_multipole(int nwords, char **words)
xyzmax = MAX(xyzmax,ypole[itype][iframe]);
if (xyzmax < 0) mpaxis[itype][iframe] = THREEFOLD;
if (mpaxis[itype][iframe] < 0) error->all(FLERR,"Mpaxis value not set");
// now reset xyz pole to positive values
if (xpole[itype][iframe] < 0) xpole[itype][iframe] = -xpole[itype][iframe];
@ -1151,7 +1151,7 @@ void PairAmoeba::file_charge_penetration(int nwords, char **words)
{
if (nwords < 4)
error->all(FLERR,"AMOEBA charge penetration line is invalid");
if (strstr(words[0],"chgpen") != words[0])
if (strstr(words[0],"chgpen") != words[0])
error->all(FLERR,"AMOEBA charge penetration line is invalid");
int iclass = utils::inumeric(FLERR,words[1],true,lmp);
@ -1172,7 +1172,7 @@ void PairAmoeba::file_dippolar(int nwords, char **words)
if (nwords < nparams)
error->all(FLERR,"AMOEBA dipole polarizability line is invalid");
if (strstr(words[0],"polarize") != words[0])
if (strstr(words[0],"polarize") != words[0])
error->all(FLERR,"AMOEBA dipole polarizability line is invalid");
int itype = utils::inumeric(FLERR,words[1],true,lmp);
@ -1193,7 +1193,7 @@ void PairAmoeba::file_dippolar(int nwords, char **words)
npolgroup[itype] = ngroup;
for (int igroup = 0; igroup < ngroup; igroup++)
polgroup[itype][igroup] =
polgroup[itype][igroup] =
utils::inumeric(FLERR,words[nparams+igroup],true,lmp);
}
@ -1203,7 +1203,7 @@ void PairAmoeba::file_charge_transfer(int nwords, char **words)
{
if (nwords < 4)
error->all(FLERR,"AMOEBA charge transfer line is invalid");
if (strstr(words[0],"chgtrn") != words[0])
if (strstr(words[0],"chgtrn") != words[0])
error->all(FLERR,"AMOEBA charge transfer line is invalid");
int iclass = utils::inumeric(FLERR,words[1],true,lmp);
@ -1223,7 +1223,7 @@ void PairAmoeba::initialize_type_class()
n_amtype = n_amclass = 0;
max_amtype = max_amclass = 0;
nvdwl_pair = max_vdwl_pair = 0;
// per type data
amtype_defined = NULL;

16
src/AMOEBA/amoeba_hal.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -92,7 +92,7 @@ void PairAmoeba::hal()
factor_hal = special_hal[special_which];
if (factor_hal == 0.0) continue;
j &= NEIGHMASK15;
xr = xi - xred[j][0];
yr = yi - xred[j][1];
zr = zi - xred[j][2];
@ -141,7 +141,7 @@ void PairAmoeba::hal()
}
ehal += e;
// find the chain rule terms for derivative components
de = de / rik;
@ -151,13 +151,13 @@ void PairAmoeba::hal()
// increment the total van der Waals energy and derivatives
// if jv < 0, trigger an error, needed H-bond partner is missing
iv = ired2local[i];
jv = ired2local[j];
if (jv < 0)
error->one(FLERR,"AMOEBA hal cannot find H bond partner - "
"ghost comm is too short");
error->one(FLERR,"AMOEBA hal cannot find H bond partner - "
"ghost comm is too short");
if (i == iv) {
f[i][0] -= dedx;
f[i][1] -= dedy;
@ -205,7 +205,7 @@ void PairAmoeba::hal()
// energy = e
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}

424
src/AMOEBA/amoeba_induce.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -69,7 +69,7 @@ void PairAmoeba::induce()
// set cutoffs, taper coeffs, and PME params
// create qfac here, free at end of polar()
if (use_ewald) {
choose(POLAR_LONG);
int nmine = p_kspace->nfft_owned;
@ -105,7 +105,7 @@ void PairAmoeba::induce()
memory->create(usump,nlocal,3,"ameoba/induce:usump");
// get the electrostatic field due to permanent multipoles
dfield0c(field,fieldp);
// reverse comm to sum field,fieldp from ghost atoms to owned atoms
@ -119,11 +119,11 @@ void PairAmoeba::induce()
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("AAA FIELD atom %d: field %g %g %g: fieldp %g %g %g\n",
atom->tag[i],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
atom->tag[i],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
*/
// set induced dipoles to polarizability times direct field
for (i = 0; i < nlocal; i++) {
@ -144,17 +144,17 @@ void PairAmoeba::induce()
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("AAA UDIR atom %d: udir %g %g %g: udirp %g %g %g\n",
atom->tag[i],
DEBYE*udir[i][0],DEBYE*udir[i][1],DEBYE*udir[i][2],
DEBYE*udirp[i][0],DEBYE*udirp[i][1],DEBYE*udirp[i][2]);
atom->tag[i],
DEBYE*udir[i][0],DEBYE*udir[i][1],DEBYE*udir[i][2],
DEBYE*udirp[i][0],DEBYE*udirp[i][1],DEBYE*udirp[i][2]);
*/
// get induced dipoles via the OPT extrapolation method
// NOTE: any way to rewrite these loops to avoid allocating
// uopt,uoptp with a optorder+1 dimension, just optorder ??
// since no need to store optorder+1 values after these loops
if (poltyp == OPT) {
if (poltyp == OPT) {
for (i = 0; i < nlocal; i++) {
for (j = 0; j < 3; j++) {
uopt[i][0][j] = udir[i][j];
@ -174,7 +174,7 @@ void PairAmoeba::induce()
comm->reverse_comm(this);
for (i = 0; i < nlocal; i++) {
itype = amtype[i];
itype = amtype[i];
for (j = 0; j < 3; j++) {
uopt[i][m][j] = polarity[itype] * field[i][j];
uoptp[i][m][j] = polarity[itype] * fieldp[i][j];
@ -263,15 +263,15 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("UFIELD atom %d: uind %g %g %g uinp %g %g %g "
"field %g %g %g: fieldp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
printf("UFIELD atom %d: uind %g %g %g uinp %g %g %g "
"field %g %g %g: fieldp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
*/
// set initial conjugate gradient residual and conjugate vector
for (i = 0; i < nlocal; i++) {
@ -295,7 +295,7 @@ void PairAmoeba::induce()
cfstyle = RSD;
comm->forward_comm(this);
uscale0b(BUILD,rsd,rsdp,zrsd,zrsdp);
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
crstyle = ZRSD;
comm->reverse_comm(this);
}
@ -312,15 +312,15 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("CONJ atom %d: rsd %g %g %g rsdp %g %g %g "
"conj %g %g %g: conjp %g %g %g\n",
atom->tag[i],
rsd[i][0],rsd[i][1],rsd[i][2],
rsdp[i][0],rsdp[i][1],rsdp[i][2],
conj[i][0],conj[i][1],conj[i][2],
conjp[i][0],conjp[i][1],conjp[i][2]);
printf("CONJ atom %d: rsd %g %g %g rsdp %g %g %g "
"conj %g %g %g: conjp %g %g %g\n",
atom->tag[i],
rsd[i][0],rsd[i][1],rsd[i][2],
rsdp[i][0],rsdp[i][1],rsdp[i][2],
conj[i][0],conj[i][1],conj[i][2],
conjp[i][0],conjp[i][1],conjp[i][2]);
*/
// conjugate gradient iteration of the mutual induced dipoles
while (!done) {
@ -349,15 +349,15 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("POST-COMM FIELD atom %d: field %g %g %g: fieldp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
if (atom->tag[i] == 1)
printf("POST-COMM FIELD atom %d: field %g %g %g: fieldp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2],
field[i][0],field[i][1],field[i][2],
fieldp[i][0],fieldp[i][1],fieldp[i][2]);
*/
for (i = 0; i < nlocal; i++) {
for (j = 0; j < 3; j++) {
uind[i][j] = vec[i][j];
@ -409,17 +409,17 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("POST-MPI UIND atom %d: uind %g %g %g: uinp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2]);
if (atom->tag[i] == 1)
printf("POST-MPI UIND atom %d: uind %g %g %g: uinp %g %g %g\n",
atom->tag[i],
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2]);
*/
if (pcgprec) {
cfstyle = RSD;
comm->forward_comm(this);
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
uscale0b(APPLY,rsd,rsdp,zrsd,zrsdp);
crstyle = ZRSD;
comm->reverse_comm(this);
}
@ -428,16 +428,16 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("POST-PRECOND atom %d: rsd %g %g %g: rsdp %g %g %g "
"zrsd %g %g %g: zrsdp %g %g %g\n",
atom->tag[i],
rsd[i][0],rsd[i][1],rsd[i][2],
rsdp[i][0],rsdp[i][1],rsdp[i][2],
zrsd[i][0],zrsd[i][1],zrsd[i][2],
zrsdp[i][0],zrsdp[i][1],zrsdp[i][2]);
if (atom->tag[i] == 1)
printf("POST-PRECOND atom %d: rsd %g %g %g: rsdp %g %g %g "
"zrsd %g %g %g: zrsdp %g %g %g\n",
atom->tag[i],
rsd[i][0],rsd[i][1],rsd[i][2],
rsdp[i][0],rsdp[i][1],rsdp[i][2],
zrsd[i][0],zrsd[i][1],zrsd[i][2],
zrsdp[i][0],zrsdp[i][1],zrsdp[i][2]);
*/
b = 0.0;
bp = 0.0;
@ -508,7 +508,7 @@ void PairAmoeba::induce()
uind[i][0],uind[i][1],uind[i][2],
uinp[i][0],uinp[i][1],uinp[i][2]);
*/
if (done) {
for (i = 0; i < nlocal; i++) {
term = pcgpeek * poli[i];
@ -523,25 +523,25 @@ void PairAmoeba::induce()
/*
for (i = 0; i < nlocal; i++)
if (atom->tag[i] == 1)
printf("POST-DONE UIND atom %d: uind %g %g %g: uinp %g %g %g\n",
atom->tag[i],
DEBYE*uind[i][0],DEBYE*uind[i][1],DEBYE*uind[i][2],
DEBYE*uinp[i][0],DEBYE*uinp[i][1],DEBYE*uinp[i][2]);
if (atom->tag[i] == 1)
printf("POST-DONE UIND atom %d: uind %g %g %g: uinp %g %g %g\n",
atom->tag[i],
DEBYE*uind[i][0],DEBYE*uind[i][1],DEBYE*uind[i][2],
DEBYE*uinp[i][0],DEBYE*uinp[i][1],DEBYE*uinp[i][2]);
*/
}
// terminate the calculation if dipoles failed to converge
// NOTE: could make this an error
if (iter >= maxiter || eps > epsold)
if (me == 0)
error->warning(FLERR,"AMOEBA induced dipoles did not converge");
error->warning(FLERR,"AMOEBA induced dipoles did not converge");
}
// DEBUG output to dump file
if (UIND_DEBUG)
if (UIND_DEBUG)
dump6(fp_uind,"id uindx uindy uindz uinpx uinpy uinpz",DEBYE,uind,uinp);
// deallocation of arrays
@ -622,7 +622,7 @@ void PairAmoeba::ulspred()
// derive normal equations corresponding to least squares fit
// NOTE: check all N vs N-1 indices in code from here down
} else if (polpred == LSQR) {
double ***udalt = fixudalt->tstore;
double ***upalt = fixupalt->tstore;
@ -722,7 +722,7 @@ void PairAmoeba::ufield0c(double **field, double **fieldp)
// get the real space portion of the mutual field
if (polar_rspace_flag) umutual2b(field,fieldp);
// add the self-energy portion of the mutual field
term = (4.0/3.0) * aewald*aewald*aewald / MY_PIS;
@ -731,13 +731,13 @@ void PairAmoeba::ufield0c(double **field, double **fieldp)
field[i][j] += term*uind[i][j];
fieldp[i][j] += term*uinp[i][j];
}
/*
// DEBUG
printf("UMUTUAL2B SELF i %d term %g aewald %g uind %g %g %g field %g %g %g\n",
atom->tag[i],term,aewald,
uind[i][0],uind[i][1],uind[i][2],field[i][0],field[i][1],field[i][2]);
atom->tag[i],term,aewald,
uind[i][0],uind[i][1],uind[i][2],field[i][0],field[i][1],field[i][2]);
*/
}
}
@ -748,7 +748,7 @@ void PairAmoeba::ufield0c(double **field, double **fieldp)
gradient induced dipole solver using a neighbor pair list
------------------------------------------------------------------------- */
void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
double **zrsd, double **zrsdp)
{
int i,j,k,m,itype,jtype,iclass,jclass,igroup,jgroup;
@ -806,7 +806,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
}
// zero zrsd,zrsdp for ghost atoms only
for (i = nlocal; i < nall; i++) {
for (j = 0; j < 3; j++) {
zrsd[i][j] = 0.0;
@ -824,7 +824,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK15;
j &= NEIGHMASK15;
m1 = pclist[0];
m2 = pclist[1];
@ -866,7 +866,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
itype = amtype[i];
iclass = amtype2class[itype];
igroup = amgroup[i];
jlist = firstneigh_precond[i];
jnum = numneigh_precond[i];
@ -929,7 +929,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
scale3 = factor_wscale * dmpik[2];
scale5 = factor_wscale * dmpik[4];
}
polik = poli * polarity[jtype];
rr3 = scale3 * polik / (r*r2);
rr5 = 3.0 * scale5 * polik / (r*r2*r2);
@ -942,16 +942,16 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
pclist[5] = rr5*zr*zr - rr3;
// DEBUG
/*
printf("PCLIST: ij %d %d: pc1-6 %g %g %g %g %g %g\n",
atom->tag[i],atom->tag[j],
pclist[0],pclist[1],pclist[2],pclist[3],pclist[4],pclist[5]);
atom->tag[i],atom->tag[j],
pclist[0],pclist[1],pclist[2],pclist[3],pclist[4],pclist[5]);
printf(" AMOEBA: scale35 %g %g pdamp %g thole %g "
"pgamma %g facUscale %g damp %g\n",
scale3,scale5,pdamp[jtype],thole[jtype],pgamma,factor_uscale,damp);
"pgamma %g facUscale %g damp %g\n",
scale3,scale5,pdamp[jtype],thole[jtype],pgamma,factor_uscale,damp);
*/
pclist += 6;
}
}
@ -1023,7 +1023,7 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
double **fuind,**fuinp;
double **fdip_phi1,**fdip_phi2,**fdip_sum_phi;
double **dipfield1,**dipfield2;
// return if the Ewald coefficient is zero
if (aewald < 1.0e-6) return;
@ -1061,7 +1061,7 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
*/
// gridpre = my portion of 4d grid in brick decomp w/ ghost values
double ****gridpre = (double ****) ic_kspace->zero();
// map 2 values to grid
@ -1070,11 +1070,11 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
// pre-convolution operations including forward FFT
// gridfft = my portion of complex 3d grid in FFT decomposition
double *gridfft = ic_kspace->pre_convolution();
// ---------------------
// convolution operation
// convolution operation
// ---------------------
nxlo = ic_kspace->nxlo_fft;
@ -1093,7 +1093,7 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
term = qfac[m++];
gridfft[n] *= term;
gridfft[n+1] *= term;
n += 2;
n += 2;
}
}
}
@ -1104,7 +1104,7 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
double ****gridpost = (double ****) ic_kspace->post_convolution();
// get potential
fphi_uind(gridpost,fdip_phi1,fdip_phi2,fdip_sum_phi);
// printf ("fdip_phi1_uind %g %g %g %g %g %g %g %g %g %g\n",
@ -1138,9 +1138,9 @@ void PairAmoeba::umutual1(double **field, double **fieldp)
for (i = 0; i < nlocal; i++) {
for (j = 0; j < 3; j++) {
dipfield1[i][j] = a[j][0]*fdip_phi1[i][1] +
dipfield1[i][j] = a[j][0]*fdip_phi1[i][1] +
a[j][1]*fdip_phi1[i][2] + a[j][2]*fdip_phi1[i][3];
dipfield2[i][j] = a[j][0]*fdip_phi2[i][1] +
dipfield2[i][j] = a[j][0]*fdip_phi2[i][1] +
a[j][1]*fdip_phi2[i][2] + a[j][2]*fdip_phi2[i][3];
}
}
@ -1205,15 +1205,15 @@ void PairAmoeba::umutual2b(double **field, double **fieldp)
fid[0] = tdipdip[0]*uindj[0] + tdipdip[1]*uindj[1] + tdipdip[2]*uindj[2];
fid[1] = tdipdip[1]*uindj[0] + tdipdip[3]*uindj[1] + tdipdip[4]*uindj[2];
fid[2] = tdipdip[2]*uindj[0] + tdipdip[4]*uindj[1] + tdipdip[5]*uindj[2];
fkd[0] = tdipdip[0]*uindi[0] + tdipdip[1]*uindi[1] + tdipdip[2]*uindi[2];
fkd[1] = tdipdip[1]*uindi[0] + tdipdip[3]*uindi[1] + tdipdip[4]*uindi[2];
fkd[2] = tdipdip[2]*uindi[0] + tdipdip[4]*uindi[1] + tdipdip[5]*uindi[2];
fip[0] = tdipdip[0]*uinpj[0] + tdipdip[1]*uinpj[1] + tdipdip[2]*uinpj[2];
fip[1] = tdipdip[1]*uinpj[0] + tdipdip[3]*uinpj[1] + tdipdip[4]*uinpj[2];
fip[2] = tdipdip[2]*uinpj[0] + tdipdip[4]*uinpj[1] + tdipdip[5]*uinpj[2];
fkp[0] = tdipdip[0]*uinpi[0] + tdipdip[1]*uinpi[1] + tdipdip[2]*uinpi[2];
fkp[1] = tdipdip[1]*uinpi[0] + tdipdip[3]*uinpi[1] + tdipdip[4]*uinpi[2];
fkp[2] = tdipdip[2]*uinpi[0] + tdipdip[4]*uinpi[1] + tdipdip[5]*uinpi[2];
@ -1222,21 +1222,21 @@ void PairAmoeba::umutual2b(double **field, double **fieldp)
/*
if (atom->tag[i] == 1 || atom->tag[j] == 1) {
printf ("TDIPDIP ij %d %d: tdd %g %g %g %g %g %g\n",
atom->tag[i],atom->tag[j],
tdipdip[0],tdipdip[1],tdipdip[2],
tdipdip[3],tdipdip[4],tdipdip[5]);
printf ("FIDFKD ij %d %d: fid %g %g %g fkd %g %g %g\n",
atom->tag[i],atom->tag[j],
fid[0],fid[1],fid[2],
fkd[0],fkd[1],fkd[2]);
}
*/
printf ("TDIPDIP ij %d %d: tdd %g %g %g %g %g %g\n",
atom->tag[i],atom->tag[j],
tdipdip[0],tdipdip[1],tdipdip[2],
tdipdip[3],tdipdip[4],tdipdip[5]);
printf ("FIDFKD ij %d %d: fid %g %g %g fkd %g %g %g\n",
atom->tag[i],atom->tag[j],
fid[0],fid[1],fid[2],
fkd[0],fkd[1],fkd[2]);
}
*/
tdipdip += 6;
// increment the field at each site due to this interaction
for (m = 0; m < 3; m++) {
field[i][m] += fid[m];
field[j][m] += fkd[m];
@ -1264,7 +1264,7 @@ void PairAmoeba::udirect1(double **field)
double volterm,denom;
double hsq,expterm;
double term,pterm;
// return if the Ewald coefficient is zero
if (aewald < 1.0e-6) return;
@ -1282,7 +1282,7 @@ void PairAmoeba::udirect1(double **field)
memory->create(fmp,nlocal,10,"ameoba/induce:fmp");
memory->create(cphi,nlocal,10,"ameoba/induce:cphi");
memory->create(fphi,nlocal,20,"ameoba/induce:fphi");
// FFT moduli pre-computations
// set igrid for each atom and its B-spline coeffs
@ -1308,14 +1308,14 @@ void PairAmoeba::udirect1(double **field)
cmp[i][8] = 2.0 * rpole[i][6];
cmp[i][9] = 2.0 * rpole[i][9];
}
// convert Cartesian multipoles to fractional coordinates
cmp_to_fmp(cmp,fmp);
// gridpre = my portion of 3d grid in brick decomp w/ ghost values
// zeroed by setup()
double ***gridpre = (double ***) i_kspace->zero();
// map multipole moments to grid
@ -1324,11 +1324,11 @@ void PairAmoeba::udirect1(double **field)
// pre-convolution operations including forward FFT
// gridfft = my 1d portion of complex 3d grid in FFT decomp
double *gridfft = i_kspace->pre_convolution();
// ---------------------
// convolution operation
// convolution operation
// ---------------------
nhalf1 = (nfft1+1) / 2;
@ -1341,28 +1341,28 @@ void PairAmoeba::udirect1(double **field)
nyhi = i_kspace->nyhi_fft;
nzlo = i_kspace->nzlo_fft;
nzhi = i_kspace->nzhi_fft;
m = n = 0;
for (k = nzlo; k <= nzhi; k++) {
for (j = nylo; j <= nyhi; j++) {
for (i = nxlo; i <= nxhi; i++) {
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
term = -pterm * hsq;
expterm = 0.0;
if (term > -50.0 && hsq != 0.0) {
denom = volterm*hsq*bsmod1[i]*bsmod2[j]*bsmod3[k];
expterm = exp(term) / denom;
}
qfac[m++] = expterm;
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
term = -pterm * hsq;
expterm = 0.0;
if (term > -50.0 && hsq != 0.0) {
denom = volterm*hsq*bsmod1[i]*bsmod2[j]*bsmod3[k];
expterm = exp(term) / denom;
}
qfac[m++] = expterm;
gridfft[n] *= expterm;
gridfft[n+1] *= expterm;
n += 2;
n += 2;
}
}
}
@ -1373,11 +1373,11 @@ void PairAmoeba::udirect1(double **field)
double ***gridpost = (double ***) i_kspace->post_convolution();
// get potential
fphi_mpole(gridpost,fphi);
// printf ("fphi %g %g %g %g %g %g %g %g %g %g %g %g\n",
// fphi[1][1],fphi[1][2],fphi[1][3],fphi[1][4],fphi[1][5],
// printf ("fphi %g %g %g %g %g %g %g %g %g %g %g %g\n",
// fphi[1][1],fphi[1][2],fphi[1][3],fphi[1][4],fphi[1][5],
// fphi[1][6],fphi[1][7],fphi[1][8],fphi[1][9],fphi[1][10],
// fphi[1][11],fphi[1][12]);
@ -1453,7 +1453,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
firstneigh = list->firstneigh;
// NOTE: doesn't this have a problem if aewald is tiny ??
aesq2 = 2.0 * aewald * aewald;
aesq2n = 0.0;
if (aewald > 0.0) aesq2n = 1.0 / (MY_PIS*aewald);
@ -1505,17 +1505,17 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
printf("Atom 4 is in group %d\n",igroup);
printf("Atoms in same group:");
for (int ig = 0; ig < atom->nlocal; ig++)
if (amgroup[ig] == igroup) printf(" %d",atom->tag[ig]);
if (amgroup[ig] == igroup) printf(" %d",atom->tag[ig]);
printf("\n");
}
*/
// evaluate all sites within the cutoff distance
for (jj = 0; jj < jnum; jj++) {
jextra = jlist[jj];
j = jextra & NEIGHMASK15;
xr = x[j][0] - x[i][0];
yr = x[j][1] - x[i][1];
zr = x[j][2] - x[i][2];
@ -1525,36 +1525,36 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
jtype = amtype[j];
jclass = amtype2class[jtype];
jgroup = amgroup[j];
if (amoeba) {
factor_wscale = special_polar_wscale[sbmask15(jextra)];
if (igroup == jgroup) {
factor_pscale = special_polar_piscale[sbmask15(jextra)];
factor_dscale = polar_dscale;
factor_uscale = polar_uscale;
} else {
factor_pscale = special_polar_pscale[sbmask15(jextra)];
factor_dscale = factor_uscale = 1.0;
}
factor_wscale = special_polar_wscale[sbmask15(jextra)];
if (igroup == jgroup) {
factor_pscale = special_polar_piscale[sbmask15(jextra)];
factor_dscale = polar_dscale;
factor_uscale = polar_uscale;
} else {
factor_pscale = special_polar_pscale[sbmask15(jextra)];
factor_dscale = factor_uscale = 1.0;
}
} else if (hippo) {
factor_wscale = special_polar_wscale[sbmask15(jextra)];
if (igroup == jgroup) {
factor_dscale = factor_pscale = special_polar_piscale[sbmask15(jextra)];
factor_uscale = polar_uscale;
} else {
factor_dscale = factor_pscale = special_polar_pscale[sbmask15(jextra)];
factor_uscale = 1.0;
}
factor_wscale = special_polar_wscale[sbmask15(jextra)];
if (igroup == jgroup) {
factor_dscale = factor_pscale = special_polar_piscale[sbmask15(jextra)];
factor_uscale = polar_uscale;
} else {
factor_dscale = factor_pscale = special_polar_pscale[sbmask15(jextra)];
factor_uscale = 1.0;
}
}
// DEBUG
/*
if (atom->tag[i] == 4 || atom->tag[j] == 4) {
printf("PAIR ij %d %d ij group %d %d wpdu scale %g %g %g %g\n",
atom->tag[i],atom->tag[j],igroup,jgroup,
factor_wscale,factor_pscale,factor_dscale,factor_uscale);
printf("PAIR ij %d %d ij group %d %d wpdu scale %g %g %g %g\n",
atom->tag[i],atom->tag[j],igroup,jgroup,
factor_wscale,factor_pscale,factor_dscale,factor_uscale);
}
*/
@ -1574,7 +1574,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
qkyy = rpole[j][8];
qkyz = rpole[j][9];
qkzz = rpole[j][12];
// intermediates involving moments and separation distance
dir = dix*xr + diy*yr + diz*zr;
@ -1587,7 +1587,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
qky = qkxy*xr + qkyy*yr + qkyz*zr;
qkz = qkxz*xr + qkyz*yr + qkzz*zr;
qkr = qkx*xr + qky*yr + qkz*zr;
// calculate the real space Ewald error function terms
ralpha = aewald * r;
@ -1599,7 +1599,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
aefac = aesq2 * aefac;
bn[m] = (bfac*bn[m-1]+aefac*exp2a) * rr2;
}
// find the field components for Thole polarization damping
if (amoeba) {
@ -1633,24 +1633,24 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
bcn[2] = bn[3] - (1.0-scalek*scale7)*rr7;
fid[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
fid[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
bcn[0]*dkx + 2.0*bcn[1]*qkx;
fid[1] = -yr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
fid[1] = -yr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
bcn[0]*dky + 2.0*bcn[1]*qky;
fid[2] = -zr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
fid[2] = -zr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
bcn[0]*dkz + 2.0*bcn[1]*qkz;
fkd[0] = xr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
fkd[0] = xr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
bcn[0]*dix - 2.0*bcn[1]*qix;
fkd[1] = yr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
fkd[1] = yr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
bcn[0]*diy - 2.0*bcn[1]*qiy;
fkd[2] = zr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
fkd[2] = zr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
bcn[0]*diz - 2.0*bcn[1]*qiz;
scalek = factor_pscale;
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
bcn[2] = bn[3] - (1.0-scalek*scale7)*rr7;
fip[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
fip[0] = -xr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
bcn[0]*dkx + 2.0*bcn[1]*qkx;
fip[1] = -yr*(bcn[0]*ck-bcn[1]*dkr+bcn[2]*qkr) -
bcn[0]*dky + 2.0*bcn[1]*qky;
@ -1662,7 +1662,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
bcn[0]*diy - 2.0*bcn[1]*qiy;
fkp[2] = zr*(bcn[0]*ci+bcn[1]*dir+bcn[2]*qir) -
bcn[0]*diz - 2.0*bcn[1]*qiz;
// find terms needed later to compute mutual polarization
if (poltyp != DIRECT) {
@ -1681,7 +1681,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
scalek = factor_uscale;
bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3;
bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5;
neighptr[n++] = j;
tdipdip[ndip++] = -bcn[0] + bcn[1]*xr*xr;
tdipdip[ndip++] = bcn[1]*xr*yr;
@ -1690,7 +1690,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
tdipdip[ndip++] = bcn[1]*yr*zr;
tdipdip[ndip++] = -bcn[0] + bcn[1]*zr*zr;
}
// find the field components for charge penetration damping
} else if (hippo) {
@ -1698,7 +1698,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
alphak = palpha[jclass];
valk = pval[j];
dampdir(r,alphai,alphak,dmpi,dmpk);
scalek = factor_dscale;
rr3i = bn[1] - (1.0-scalek*dmpi[2])*rr3;
rr5i = bn[2] - (1.0-scalek*dmpi[4])*rr5;
@ -1741,7 +1741,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
rr3i*diy - 2.0*rr5i*qiy;
fkp[2] = zr*(rr3*corei + rr3i*vali + rr5i*dir + rr7i*qir) -
rr3i*diz - 2.0*rr5i*qiz;
// find terms needed later to compute mutual polarization
if (poltyp != DIRECT) {
@ -1751,16 +1751,16 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
rr3ik = bn[1] - (1.0-scalek*dmpik[2])*rr3;
rr5ik = bn[2] - (1.0-scalek*dmpik[4])*rr5;
// DEBUG
// DEBUG
/*
if (atom->tag[i] == 1 || atom->tag[j] == 1)
printf ("DAMPMUT ij %d %d: bn12 %g %g dmpik-01234 %g %g %g %g %g\n",
atom->tag[i],atom->tag[j],
bn[1],bn[2],
dmpik[0],dmpik[1],dmpik[2],dmpik[3],dmpik[4]);
*/
/*
if (atom->tag[i] == 1 || atom->tag[j] == 1)
printf ("DAMPMUT ij %d %d: bn12 %g %g dmpik-01234 %g %g %g %g %g\n",
atom->tag[i],atom->tag[j],
bn[1],bn[2],
dmpik[0],dmpik[1],dmpik[2],dmpik[3],dmpik[4]);
*/
neighptr[n++] = j;
tdipdip[ndip++] = -rr3ik + rr5ik*xr*xr;
tdipdip[ndip++] = rr5ik*xr*yr;
@ -1770,9 +1770,9 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
tdipdip[ndip++] = -rr3ik + rr5ik*zr*zr;
}
}
// increment the field at each site due to this interaction
for (m = 0; m < 3; m++) {
field[i][m] += fid[m];
field[j][m] += fkd[m];
@ -1823,9 +1823,9 @@ void PairAmoeba::dampmut(double r, double alphai, double alphak, double *dmpik)
if (diff < eps) {
dampi4 = dampi2 * dampi2;
dampi5 = dampi2 * dampi3;
dmpik[2] = 1.0 - (1.0 + dampi + 0.5*dampi2 +
dmpik[2] = 1.0 - (1.0 + dampi + 0.5*dampi2 +
7.0*dampi3/48.0 + dampi4/48.0)*expi;
dmpik[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dmpik[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/144.0)*expi;
} else {
dampk2 = dampk * dampk;
@ -1836,12 +1836,12 @@ void PairAmoeba::dampmut(double r, double alphai, double alphak, double *dmpik)
termk = alphai2 / (alphai2-alphak2);
termi2 = termi * termi;
termk2 = termk * termk;
dmpik[2] = 1.0 - termi2*(1.0+dampi+0.5*dampi2)*expi -
termk2*(1.0+dampk+0.5*dampk2)*expk -
dmpik[2] = 1.0 - termi2*(1.0+dampi+0.5*dampi2)*expi -
termk2*(1.0+dampk+0.5*dampk2)*expk -
2.0*termi2*termk*(1.0+dampi)*expi - 2.0*termk2*termi*(1.0+dampk)*expk;
dmpik[4] = 1.0 - termi2*(1.0+dampi+0.5*dampi2 + dampi3/6.0)*expi -
termk2*(1.0+dampk+0.5*dampk2 + dampk3/6.00)*expk -
2.0*termi2*termk *(1.0+dampi+dampi2/3.0)*expi -
dmpik[4] = 1.0 - termi2*(1.0+dampi+0.5*dampi2 + dampi3/6.0)*expi -
termk2*(1.0+dampk+0.5*dampk2 + dampk3/6.00)*expk -
2.0*termi2*termk *(1.0+dampi+dampi2/3.0)*expi -
2.0*termk2*termi *(1.0+dampk+dampk2/3.0)*expk;
}
}
@ -1852,7 +1852,7 @@ void PairAmoeba::dampmut(double r, double alphai, double alphak, double *dmpik)
function for powers of the interatomic distance
------------------------------------------------------------------------- */
void PairAmoeba::dampdir(double r, double alphai, double alphak,
void PairAmoeba::dampdir(double r, double alphai, double alphak,
double *dmpi, double *dmpk)
{
double eps,diff;
@ -1917,7 +1917,7 @@ void PairAmoeba::cholesky(int nvar, double *a, double *b)
a--;
b--;
// Cholesky factorization to reduce A to (L)(D)(L transpose)
// L has a unit diagonal; store 1.0/D on the diagonal of A
@ -1927,19 +1927,19 @@ void PairAmoeba::cholesky(int nvar, double *a, double *b)
if (i != 1) {
ij = i;
for (j = 1; j <= im; j++) {
r = a[ij];
if (j != 1) {
ik = i;
jk = j;
jm = j - 1;
for (k = 1; k <= jm; k++) {
r = r - a[ik]*a[jk];
ik = nvar - k + ik;
jk = nvar - k + jk;
}
}
a[ij] = r;
ij = nvar - j + ij;
r = a[ij];
if (j != 1) {
ik = i;
jk = j;
jm = j - 1;
for (k = 1; k <= jm; k++) {
r = r - a[ik]*a[jk];
ik = nvar - k + ik;
jk = nvar - k + jk;
}
}
a[ij] = r;
ij = nvar - j + ij;
}
}
@ -1948,12 +1948,12 @@ void PairAmoeba::cholesky(int nvar, double *a, double *b)
kk = 1;
ik = i;
for (k = 1; k <= im; k++) {
s = a[ik];
t = s * a[kk];
a[ik] = t;
r = r - s*t;
ik = nvar - k + ik;
kk = nvar - k + 1 + kk;
s = a[ik];
t = s * a[kk];
a[ik] = t;
r = r - s*t;
ik = nvar - k + ik;
kk = nvar - k + 1 + kk;
}
}
@ -1969,8 +1969,8 @@ void PairAmoeba::cholesky(int nvar, double *a, double *b)
im = i - 1;
r = b[i];
for (k = 1; k <= im; k++) {
r = r - b[k]*a[ik];
ik = nvar - k + ik;
r = r - b[k]*a[ik];
ik = nvar - k + ik;
}
b[i] = r;
}
@ -1986,8 +1986,8 @@ void PairAmoeba::cholesky(int nvar, double *a, double *b)
im = i + 1;
ki = ii + 1;
for (k = im; k <= nvar; k++) {
r = r - a[ki]*b[k];
ki = ki + 1;
r = r - a[ki]*b[k];
ki = ki + 1;
}
}
b[i] = r;

372
src/AMOEBA/amoeba_kspace.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -61,7 +61,7 @@ void PairAmoeba::moduli()
int maxfft = MAX(nfft1,nfft2);
maxfft = MAX(maxfft,nfft3);
double *array = new double[bsorder];
double *bsarray = new double[maxfft];
@ -79,7 +79,7 @@ void PairAmoeba::moduli()
dftmod(bsmod3,bsarray,nfft3,bsorder);
// perform deallocation of local arrays
delete [] array;
delete [] bsarray;
}
@ -104,7 +104,7 @@ void PairAmoeba::bspline(double x, int n, double *c)
/*for (k = 2; k < n; k++) {
denom = 1.0 / (k-1);
c[k] = x * c[k-1] * denom;
for (i = 0; i < k-1; i++)
for (i = 0; i < k-1; i++)
c[k-i] = ((x+i)*c[k-i-1] + (k-i-x)*c[k-i]) * denom;
c[0] = (1.0-x) * c[0] * denom;
}*/
@ -148,15 +148,15 @@ void PairAmoeba::dftmod(double *bsmod, double *bsarray, int nfft, int order)
bsmod[i] = sum1*sum1 + sum2*sum2;
//printf("BSMOD AAA i %d bsmod %g\n",i,bsmod[i]);
}
// fix for exponential Euler spline interpolation failure
eps = 1.0e-7;
if (bsmod[0] < eps) bsmod[0] = 0.5 * bsmod[1];
for (i = 1; i < nfft-1; i++)
for (i = 1; i < nfft-1; i++)
if (bsmod[i] < eps) bsmod[i] = 0.5 * (bsmod[i-1] + bsmod[i+1]);
if (bsmod[nfft-1] < eps) bsmod[nfft-1] = 0.5 * bsmod[nfft-2];
// compute and apply the optimal zeta coefficient
jcut = 50;
@ -170,14 +170,14 @@ void PairAmoeba::dftmod(double *bsmod, double *bsarray, int nfft, int order)
sum2 = 1.0;
factor = MY_PI * k / nfft;
for (j = 1; j <= jcut; j++) {
arg = factor / (factor + MY_PI*j);
sum1 += pow(arg,order);
sum2 += pow(arg,order2);
arg = factor / (factor + MY_PI*j);
sum1 += pow(arg,order);
sum2 += pow(arg,order2);
}
for (j = 1; j <= jcut; j++) {
arg = factor / (factor - MY_PI*j);
sum1 += pow(arg,order);
sum2 += pow(arg,order2);
arg = factor / (factor - MY_PI*j);
sum1 += pow(arg,order);
sum2 += pow(arg,order2);
}
zeta = sum2 / sum1;
}
@ -197,7 +197,7 @@ void PairAmoeba::bspline_fill()
double xi,yi,zi;
double w,fr,eps;
double lamda[3];
int nlocal = atom->nlocal;
double **x = atom->x;
@ -258,7 +258,7 @@ void PairAmoeba::bsplgen(double w, double **thetai)
double denom;
level = 4;
// initialization to get to 2nd order recursion
bsbuild[1][1] = w;
@ -278,7 +278,7 @@ void PairAmoeba::bsplgen(double w, double **thetai)
bsbuild[i-1][i-1] = denom * w * bsbuild[k-1][k-1];
for (j = 1; j <= i-2; j++) {
bsbuild[i-j-1][i-1] = denom *
((w+j)*bsbuild[i-j-2][k-1] + (i-j-w)*bsbuild[i-j-1][k-1]);
((w+j)*bsbuild[i-j-2][k-1] + (i-j-w)*bsbuild[i-j-1][k-1]);
}
bsbuild[0][i-1] = denom * (1.0-w) * bsbuild[0][k-1];
}
@ -290,7 +290,7 @@ void PairAmoeba::bsplgen(double w, double **thetai)
for (int i = bsorder-1; i >= 2; i--)
bsbuild[i-1][k-1] = bsbuild[i-2][k-1] - bsbuild[i-1][k-1];
bsbuild[0][k-1] = -bsbuild[0][k-1];
// get coefficients for the B-spline second derivative
if (level == 4) {
@ -351,7 +351,7 @@ void PairAmoeba::cmp_to_fmp(double **cmp, double **fmp)
for (j = 1; j < 4; j++) {
fmp[i][j] = 0.0;
for (k = 1; k < 4; k++)
fmp[i][j] += ctf[j][k] * cmp[i][k];
fmp[i][j] += ctf[j][k] * cmp[i][k];
}
for (j = 4; j < 10; j++) {
fmp[i][j] = 0.0;
@ -439,12 +439,12 @@ void PairAmoeba::fphi_to_cphi(double **fphi, double **cphi)
for (j = 1; j < 4; j++) {
cphi[i][j] = 0.0;
for (k = 1; k < 4; k++)
cphi[i][j] += ftc[j][k] * fphi[i][k];
cphi[i][j] += ftc[j][k] * fphi[i][k];
}
for (j = 4; j < 10; j++) {
cphi[i][j] = 0.0;
for (k = 4; k < 10; k++)
cphi[i][j] += ftc[j][k] * fphi[i][k];
cphi[i][j] += ftc[j][k] * fphi[i][k];
}
}
}
@ -530,9 +530,9 @@ void PairAmoeba::grid_mpole(double **fmp, double ***grid)
// spread the permanent multipole moments onto the grid
int nlocal = atom->nlocal;
for (m = 0; m < nlocal; m++) {
k = igrid[m][2] - nlpts;
for (kb = 0; kb < bsorder; kb++) {
v0 = thetai3[m][kb][0];
@ -541,43 +541,43 @@ void PairAmoeba::grid_mpole(double **fmp, double ***grid)
j = igrid[m][1] - nlpts;
for (jb = 0; jb < bsorder; jb++) {
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
term0 = fmp[m][0]*u0*v0 + fmp[m][2]*u1*v0 + fmp[m][3]*u0*v1 +
fmp[m][5]*u2*v0 + fmp[m][6]*u0*v2 + fmp[m][9]*u1*v1;
term1 = fmp[m][1]*u0*v0 + fmp[m][7]*u1*v0 + fmp[m][8]*u0*v1;
term2 = fmp[m][4]*u0*v0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
t1 = thetai1[m][ib][1];
t2 = thetai1[m][ib][2];
grid[k][j][i] += term0*t0 + term1*t1 + term2*t2;
// if (m == 0) {
// int istencil = kb*bsorder*bsorder + jb*bsorder + ib + 1;
// printf("GRIDMPOLE iStencil %d atomID %d igrid %d %d %d "
// "ibjbkb %d %d %d ijk %d %d %d "
// "th1 %g %g %g th2 %g %g %g th3 %g %g %g "
// "fmp0-9 %e %e %e %e %e %e %e %e %e %e "
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
term0 = fmp[m][0]*u0*v0 + fmp[m][2]*u1*v0 + fmp[m][3]*u0*v1 +
fmp[m][5]*u2*v0 + fmp[m][6]*u0*v2 + fmp[m][9]*u1*v1;
term1 = fmp[m][1]*u0*v0 + fmp[m][7]*u1*v0 + fmp[m][8]*u0*v1;
term2 = fmp[m][4]*u0*v0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
t1 = thetai1[m][ib][1];
t2 = thetai1[m][ib][2];
grid[k][j][i] += term0*t0 + term1*t1 + term2*t2;
// if (m == 0) {
// int istencil = kb*bsorder*bsorder + jb*bsorder + ib + 1;
// printf("GRIDMPOLE iStencil %d atomID %d igrid %d %d %d "
// "ibjbkb %d %d %d ijk %d %d %d "
// "th1 %g %g %g th2 %g %g %g th3 %g %g %g "
// "fmp0-9 %e %e %e %e %e %e %e %e %e %e "
// "term012 %e %e %e "
// "gridvalue %g\n",
// istencil,atom->tag[m],
// igrid[m][0],igrid[m][1],igrid[m][2],
// ib,jb,kb,i,j,k,
// t0,t1,t2,u0,u1,u2,v0,v1,v2,
// fmp[m][0],fmp[m][1],fmp[m][2],
// fmp[m][3],fmp[m][4],fmp[m][5],
// fmp[m][6],fmp[m][7],fmp[m][8],fmp[m][9],
// "gridvalue %g\n",
// istencil,atom->tag[m],
// igrid[m][0],igrid[m][1],igrid[m][2],
// ib,jb,kb,i,j,k,
// t0,t1,t2,u0,u1,u2,v0,v1,v2,
// fmp[m][0],fmp[m][1],fmp[m][2],
// fmp[m][3],fmp[m][4],fmp[m][5],
// fmp[m][6],fmp[m][7],fmp[m][8],fmp[m][9],
// term0,term1,term2,
// term0*t0 + term1*t1 + term2*t2);
// }
i++;
}
j++;
// term0*t0 + term1*t1 + term2*t2);
// }
i++;
}
j++;
}
k++;
}
@ -608,9 +608,9 @@ void PairAmoeba::fphi_mpole(double ***grid, double **fphi)
int nrpts = bsorder - nlpts - 1;
// extract the permanent multipole field at each site
int nlocal = atom->nlocal;
for (m = 0; m < nlocal; m++) {
tuv000 = 0.0;
tuv001 = 0.0;
@ -649,39 +649,39 @@ void PairAmoeba::fphi_mpole(double ***grid, double **fphi)
tu21 = 0.0;
tu12 = 0.0;
tu03 = 0.0;
j = igrid[m][1] - nlpts;
for (jb = 0; jb < bsorder; jb++) {
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
u3 = thetai2[m][jb][3];
t0 = 0.0;
t1 = 0.0;
t2 = 0.0;
t3 = 0.0;
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
u3 = thetai2[m][jb][3];
t0 = 0.0;
t1 = 0.0;
t2 = 0.0;
t3 = 0.0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
tq = grid[k][j][i];
t0 += tq*thetai1[m][ib][0];
t1 += tq*thetai1[m][ib][1];
t2 += tq*thetai1[m][ib][2];
t3 += tq*thetai1[m][ib][3];
i++;
}
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
tq = grid[k][j][i];
t0 += tq*thetai1[m][ib][0];
t1 += tq*thetai1[m][ib][1];
t2 += tq*thetai1[m][ib][2];
t3 += tq*thetai1[m][ib][3];
i++;
}
tu00 += t0*u0;
tu10 += t1*u0;
tu01 += t0*u1;
tu20 += t2*u0;
tu11 += t1*u1;
tu02 += t0*u2;
tu30 += t3*u0;
tu21 += t2*u1;
tu12 += t1*u2;
tu03 += t0*u3;
j++;
tu00 += t0*u0;
tu10 += t1*u0;
tu01 += t0*u1;
tu20 += t2*u0;
tu11 += t1*u1;
tu02 += t0*u2;
tu30 += t3*u0;
tu21 += t2*u1;
tu12 += t1*u2;
tu03 += t0*u3;
j++;
}
tuv000 += tu00*v0;
@ -749,7 +749,7 @@ void PairAmoeba::grid_uind(double **fuind, double **fuinp, double ****grid)
// put the induced dipole moments onto the grid
int nlocal = atom->nlocal;
for (m = 0; m < nlocal; m++) {
k = igrid[m][2] - nlpts;
@ -759,23 +759,23 @@ void PairAmoeba::grid_uind(double **fuind, double **fuinp, double ****grid)
j = igrid[m][1] - nlpts;
for (jb = 0; jb < bsorder; jb++) {
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
term01 = fuind[m][1]*u1*v0 + fuind[m][2]*u0*v1;
term11 = fuind[m][0]*u0*v0;
term02 = fuinp[m][1]*u1*v0 + fuinp[m][2]*u0*v1;
term12 = fuinp[m][0]*u0*v0;
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
term01 = fuind[m][1]*u1*v0 + fuind[m][2]*u0*v1;
term11 = fuind[m][0]*u0*v0;
term02 = fuinp[m][1]*u1*v0 + fuinp[m][2]*u0*v1;
term12 = fuinp[m][0]*u0*v0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
t1 = thetai1[m][ib][1];
grid[k][j][i][0] += term01*t0 + term11*t1;
grid[k][j][i][1] += term02*t0 + term12*t1;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
t1 = thetai1[m][ib][1];
grid[k][j][i][0] += term01*t0 + term11*t1;
grid[k][j][i][1] += term02*t0 + term12*t1;
//printf("gridcheck %g %g \n",grid[k][j][i][0],grid[k][j][i][0]);
i++;
}
j++;
i++;
}
j++;
}
k++;
}
@ -783,12 +783,12 @@ void PairAmoeba::grid_uind(double **fuind, double **fuinp, double ****grid)
}
/* ----------------------------------------------------------------------
fphi_uind = induced potential from grid
fphi_uind = induced potential from grid
fphi_uind extracts the induced dipole potential from the particle mesh Ewald grid
------------------------------------------------------------------------- */
void PairAmoeba::fphi_uind(double ****grid, double **fdip_phi1,
double **fdip_phi2, double **fdip_sum_phi)
double **fdip_phi2, double **fdip_sum_phi)
{
int i,j,k,m,ib,jb,kb;
double v0,v1,v2,v3;
@ -820,7 +820,7 @@ void PairAmoeba::fphi_uind(double ****grid, double **fdip_phi1,
// extract the permanent multipole field at each site
int nlocal = atom->nlocal;
for (m = 0; m < nlocal; m++) {
tuv100_1 = 0.0;
tuv010_1 = 0.0;
@ -892,60 +892,60 @@ void PairAmoeba::fphi_uind(double ****grid, double **fdip_phi1,
j = igrid[m][1] - nlpts;
for (jb = 0; jb < bsorder; jb++) {
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
u3 = thetai2[m][jb][3];
t0_1 = 0.0;
t1_1 = 0.0;
t2_1 = 0.0;
t0_2 = 0.0;
t1_2 = 0.0;
t2_2 = 0.0;
t3 = 0.0;
u0 = thetai2[m][jb][0];
u1 = thetai2[m][jb][1];
u2 = thetai2[m][jb][2];
u3 = thetai2[m][jb][3];
t0_1 = 0.0;
t1_1 = 0.0;
t2_1 = 0.0;
t0_2 = 0.0;
t1_2 = 0.0;
t2_2 = 0.0;
t3 = 0.0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
tq_1 = grid[k][j][i][0];
tq_2 = grid[k][j][i][1];
t0_1 += tq_1*thetai1[m][ib][0];
t1_1 += tq_1*thetai1[m][ib][1];
t2_1 += tq_1*thetai1[m][ib][2];
t0_2 += tq_2*thetai1[m][ib][0];
t1_2 += tq_2*thetai1[m][ib][1];
t2_2 += tq_2*thetai1[m][ib][2];
t3 += (tq_1+tq_2)*thetai1[m][ib][3];
i++;
}
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
tq_1 = grid[k][j][i][0];
tq_2 = grid[k][j][i][1];
t0_1 += tq_1*thetai1[m][ib][0];
t1_1 += tq_1*thetai1[m][ib][1];
t2_1 += tq_1*thetai1[m][ib][2];
t0_2 += tq_2*thetai1[m][ib][0];
t1_2 += tq_2*thetai1[m][ib][1];
t2_2 += tq_2*thetai1[m][ib][2];
t3 += (tq_1+tq_2)*thetai1[m][ib][3];
i++;
}
tu00_1 += t0_1*u0;
tu10_1 += t1_1*u0;
tu01_1 += t0_1*u1;
tu20_1 += t2_1*u0;
tu11_1 += t1_1*u1;
tu02_1 += t0_1*u2;
tu00_2 += t0_2*u0;
tu10_2 += t1_2*u0;
tu01_2 += t0_2*u1;
tu20_2 += t2_2*u0;
tu11_2 += t1_2*u1;
tu02_2 += t0_2*u2;
t0 = t0_1 + t0_2;
t1 = t1_1 + t1_2;
t2 = t2_1 + t2_2;
tu00 += t0*u0;
tu10 += t1*u0;
tu01 += t0*u1;
tu20 += t2*u0;
tu11 += t1*u1;
tu02 += t0*u2;
tu30 += t3*u0;
tu21 += t2*u1;
tu12 += t1*u2;
tu03 += t0*u3;
j++;
tu00_1 += t0_1*u0;
tu10_1 += t1_1*u0;
tu01_1 += t0_1*u1;
tu20_1 += t2_1*u0;
tu11_1 += t1_1*u1;
tu02_1 += t0_1*u2;
tu00_2 += t0_2*u0;
tu10_2 += t1_2*u0;
tu01_2 += t0_2*u1;
tu20_2 += t2_2*u0;
tu11_2 += t1_2*u1;
tu02_2 += t0_2*u2;
t0 = t0_1 + t0_2;
t1 = t1_1 + t1_2;
t2 = t2_1 + t2_2;
tu00 += t0*u0;
tu10 += t1*u0;
tu01 += t0*u1;
tu20 += t2*u0;
tu11 += t1*u1;
tu02 += t0*u2;
tu30 += t3*u0;
tu21 += t2*u1;
tu12 += t1*u2;
tu03 += t0*u3;
j++;
}
tuv100_1 += tu10_1*v0;
tuv010_1 += tu01_1*v0;
tuv001_1 += tu00_1*v1;
@ -997,7 +997,7 @@ void PairAmoeba::fphi_uind(double ****grid, double **fdip_phi1,
fdip_phi1[m][7] = tuv110_1;
fdip_phi1[m][8] = tuv101_1;
fdip_phi1[m][9] = tuv011_1;
fdip_phi2[m][0] = 0.0;
fdip_phi2[m][1] = tuv100_2;
fdip_phi2[m][2] = tuv010_2;
@ -1008,7 +1008,7 @@ void PairAmoeba::fphi_uind(double ****grid, double **fdip_phi1,
fdip_phi2[m][7] = tuv110_2;
fdip_phi2[m][8] = tuv101_2;
fdip_phi2[m][9] = tuv011_2;
fdip_sum_phi[m][0] = tuv000;
fdip_sum_phi[m][1] = tuv100;
fdip_sum_phi[m][2] = tuv010;
@ -1049,7 +1049,7 @@ void PairAmoeba::grid_disp(double ***grid)
// put the dispersion sites onto the grid
int nlocal = atom->nlocal;
for (m = 0; m < nlocal; m++) {
itype = amtype[m];
iclass = amtype2class[itype];
@ -1057,19 +1057,19 @@ void PairAmoeba::grid_disp(double ***grid)
k = igrid[m][2] - nlpts;
for (kb = 0; kb < bsorder; kb++) {
v0 = thetai3[m][kb][0] * csix[iclass];
j = igrid[m][1] - nlpts;
for (jb = 0; jb < bsorder; jb++) {
u0 = thetai2[m][jb][0];
term = v0 * u0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
grid[k][j][i] += term*t0;
i++;
}
j++;
u0 = thetai2[m][jb][0];
term = v0 * u0;
i = igrid[m][0] - nlpts;
for (ib = 0; ib < bsorder; ib++) {
t0 = thetai1[m][ib][0];
grid[k][j][i] += term*t0;
i++;
}
j++;
}
k++;
}
@ -1091,19 +1091,19 @@ void PairAmoeba::kewald()
double size,slope;
// use_ewald is for multipole and induce and polar
if (!use_ewald) aeewald = apewald = 0.0;
if (use_ewald) {
if (!aeewald_key) aeewald = ewaldcof(ewaldcut);
if (!apewald_key) {
apewald = aeewald;
size = MIN(domain->xprd,domain->yprd);
size = MIN(size,domain->zprd);
if (size < 6.0) {
slope = (1.0-apewald) / 2.0;
apewald += slope*(6.0-size);
slope = (1.0-apewald) / 2.0;
apewald += slope*(6.0-size);
}
}
@ -1122,14 +1122,14 @@ void PairAmoeba::kewald()
while (!factorable(nefft2)) nefft2++;
while (!factorable(nefft3)) nefft3++;
}
// use_dewald is for dispersion
if (!use_dewald) adewald = 0.0;
if (use_dewald) {
if (!adewald_key) adewald = ewaldcof(dispcut);
if (!dpmegrid_key) {
delta = 1.0e-8;
ddens = 0.8;
@ -1140,7 +1140,7 @@ void PairAmoeba::kewald()
// increase grid sizes until factorable by 2,3,5
// NOTE: also worry about satisfying Tinker minfft ?
while (!factorable(ndfft1)) ndfft1++;
while (!factorable(ndfft2)) ndfft3++;
while (!factorable(ndfft3)) ndfft3++;
@ -1160,12 +1160,12 @@ void PairAmoeba::kewald()
if (use_dewald) nfft1 = MAX(nfft1,ndfft1);
if (use_dewald) nfft2 = MAX(nfft2,ndfft2);
if (use_dewald) nfft3 = MAX(nfft3,ndfft3);
bsordermax = 0;
if (use_ewald) bsordermax = bseorder;
if (use_ewald) bsordermax = MAX(bsordermax,bsporder);
if (use_dewald) bsordermax = MAX(bsordermax,bsdorder);
memory->create(bsmod1,nfft1,"amoeba:bsmod1");
memory->create(bsmod2,nfft2,"amoeba:bsmod2");
memory->create(bsmod3,nfft3,"amoeba:bsmod3");
@ -1203,7 +1203,7 @@ double PairAmoeba::ewaldcof(double cutoff)
}
// use a binary search to refine the coefficient
k = i + 60;
xlo = 0.0;
xhi = x;
@ -1214,7 +1214,7 @@ double PairAmoeba::ewaldcof(double cutoff)
if (ratio >= eps) xlo = x;
else xhi = x;
}
return x;
}

258
src/AMOEBA/amoeba_multipole.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -111,7 +111,7 @@ void PairAmoeba::multipole()
qizz = rpole[i][12];
cii = ci*ci;
dii = dix*dix + diy*diy + diz*diz;
qii = 2.0*(qixy*qixy+qixz*qixz+qiyz*qiyz) +
qii = 2.0*(qixy*qixy+qixz*qixz+qiyz*qiyz) +
qixx*qixx + qiyy*qiyy + qizz*qizz;
e = fterm * (cii + term*(dii/3.0+2.0*term*qii/5.0));
empole += e;
@ -214,7 +214,7 @@ void PairAmoeba::multipole_real()
//int count = 0;
//int imin,imax;
// compute the real space portion of the Ewald summation
for (ii = 0; ii < inum; ii++) {
@ -257,7 +257,7 @@ void PairAmoeba::multipole_real()
if (r2 > off2) continue;
// DEBUG
//imin = MIN(atom->tag[i],atom->tag[j]);
//imax = MAX(atom->tag[i],atom->tag[j]);
@ -292,7 +292,7 @@ void PairAmoeba::multipole_real()
qik = qix*qkx + qiy*qky + qiz*qkz;
diqk = dix*qkx + diy*qky + diz*qkz;
dkqi = dkx*qix + dky*qiy + dkz*qiz;
qiqk = 2.0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) +
qiqk = 2.0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) +
qixx*qkxx + qiyy*qkyy + qizz*qkzz;
// additional intermediates involving moments and distance
@ -339,11 +339,11 @@ void PairAmoeba::multipole_real()
dkqirx = dkqiz*yr - dkqiy*zr;
dkqiry = dkqix*zr - dkqiz*xr;
dkqirz = dkqiy*xr - dkqix*yr;
dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy -
dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy -
2.0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz - qixz*qkxy-qiyz*qkyy-qizz*qkyz);
dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz -
dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz -
2.0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz - qixx*qkxz-qixy*qkyz-qixz*qkzz);
dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix -
dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix -
2.0*(qixx*qkxy+qixy*qkyy+qixz*qkyz - qixy*qkxx-qiyy*qkxy-qiyz*qkxz);
// get reciprocal distance terms for this interaction
@ -377,11 +377,11 @@ void PairAmoeba::multipole_real()
alphak = palpha[jclass];
valk = pval[j];
/*
printf("HIPPO MPOLE ij %d %d: pcore/alpha/val I %g %g %g: J %g %g %g\n",
atom->tag[i],atom->tag[j],corei,alphai,vali,corek,alphak,valk);
*/
/*
printf("HIPPO MPOLE ij %d %d: pcore/alpha/val I %g %g %g: J %g %g %g\n",
atom->tag[i],atom->tag[j],corei,alphai,vali,corek,alphak,valk);
*/
term1 = corei*corek;
term1i = corek*vali;
term2i = corek*dir;
@ -412,16 +412,16 @@ void PairAmoeba::multipole_real()
rr11ik = bn[5] - (1.0-scalek*dmpij[10])*rr11;
rr1 = bn[0] - (1.0-scalek)*rr1;
rr3 = bn[1] - (1.0-scalek)*rr3;
e = term1*rr1 + term4ik*rr7ik + term5ik*rr9ik +
term1i*rr1i + term1k*rr1k + term1ik*rr1ik +
term2i*rr3i + term2k*rr3k + term2ik*rr3ik +
e = term1*rr1 + term4ik*rr7ik + term5ik*rr9ik +
term1i*rr1i + term1k*rr1k + term1ik*rr1ik +
term2i*rr3i + term2k*rr3k + term2ik*rr3ik +
term3i*rr5i + term3k*rr5k + term3ik*rr5ik;
// find damped multipole intermediates for force and torque
de = term1*rr3 + term4ik*rr9ik + term5ik*rr11ik +
term1i*rr3i + term1k*rr3k + term1ik*rr3ik +
term2i*rr5i + term2k*rr5k + term2ik*rr5ik +
de = term1*rr3 + term4ik*rr9ik + term5ik*rr11ik +
term1i*rr3i + term1k*rr3k + term1ik*rr3ik +
term2i*rr5i + term2k*rr5k + term2ik*rr5ik +
term3i*rr7i + term3k*rr7k + term3ik*rr7ik;
term1 = -corek*rr3i - valk*rr3ik + dkr*rr5ik - qkr*rr7ik;
term2 = corei*rr3k + vali*rr3ik + dir*rr5ik + qir*rr7ik;
@ -467,44 +467,44 @@ void PairAmoeba::multipole_real()
count++;
if (imin == 68 && imax == 1021) {
//printf("AAA %d %d %16.12g\n",imin,imax,e);
printf("AAA %d: %d %d: %d %d: %d: %16.12g\n",
me,atom->tag[i],atom->tag[j],i,j,atom->nlocal,e);
printf("XYZ %g %g %g: %g %g %g\n",
x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2]);
//printf("AAA %d %d %16.12g\n",imin,imax,e);
printf("AAA %d: %d %d: %d %d: %d: %16.12g\n",
me,atom->tag[i],atom->tag[j],i,j,atom->nlocal,e);
printf("XYZ %g %g %g: %g %g %g\n",
x[i][0],x[i][1],x[i][2],x[j][0],x[j][1],x[j][2]);
}
*/
/*
if (atom->tag[i] == 1 || atom->tag[j] == 1) {
printf("MPOLE %d %d %d: %15.12g %15.12g %g\n",
atom->tag[i],atom->tag[j],j,r,e,factor_mpole);
printf(" BN: %g %g %g: %g %g %g\n",bn[0],bn[1],bn[2],bn[3],bn[4],bn[5]);
printf("MPOLE %d %d %d: %15.12g %15.12g %g\n",
atom->tag[i],atom->tag[j],j,r,e,factor_mpole);
printf(" BN: %g %g %g: %g %g %g\n",bn[0],bn[1],bn[2],bn[3],bn[4],bn[5]);
}
*/
// compute the force components for this interaction
frcx = de*xr + term1*dix + term2*dkx + term3*(diqkx-dkqix) +
frcx = de*xr + term1*dix + term2*dkx + term3*(diqkx-dkqix) +
term4*qix + term5*qkx + term6*(qixk+qkxi);
frcy = de*yr + term1*diy + term2*dky + term3*(diqky-dkqiy) +
frcy = de*yr + term1*diy + term2*dky + term3*(diqky-dkqiy) +
term4*qiy + term5*qky + term6*(qiyk+qkyi);
frcz = de*zr + term1*diz + term2*dkz + term3*(diqkz-dkqiz) +
frcz = de*zr + term1*diz + term2*dkz + term3*(diqkz-dkqiz) +
term4*qiz + term5*qkz + term6*(qizk+qkzi);
// compute the torque components for this interaction
ttmi[0] = -rr3*dikx + term1*dirx + term3*(dqikx+dkqirx) -
ttmi[0] = -rr3*dikx + term1*dirx + term3*(dqikx+dkqirx) -
term4*qirx - term6*(qikrx+qikx);
ttmi[1] = -rr3*diky + term1*diry + term3*(dqiky+dkqiry) -
ttmi[1] = -rr3*diky + term1*diry + term3*(dqiky+dkqiry) -
term4*qiry - term6*(qikry+qiky);
ttmi[2] = -rr3*dikz + term1*dirz + term3*(dqikz+dkqirz) -
ttmi[2] = -rr3*dikz + term1*dirz + term3*(dqikz+dkqirz) -
term4*qirz - term6*(qikrz+qikz);
ttmk[0] = rr3*dikx + term2*dkrx - term3*(dqikx+diqkrx) -
ttmk[0] = rr3*dikx + term2*dkrx - term3*(dqikx+diqkrx) -
term5*qkrx - term6*(qkirx-qikx);
ttmk[1] = rr3*diky + term2*dkry - term3*(dqiky+diqkry) -
ttmk[1] = rr3*diky + term2*dkry - term3*(dqiky+diqkry) -
term5*qkry - term6*(qkiry-qiky);
ttmk[2] = rr3*dikz + term2*dkrz - term3*(dqikz+diqkrz) -
ttmk[2] = rr3*dikz + term2*dkrz - term3*(dqikz+diqkrz) -
term5*qkrz - term6*(qkirz-qikz);
// increment force-based gradient and torque on first site
@ -544,7 +544,7 @@ void PairAmoeba::multipole_real()
// energy = e
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}
@ -558,7 +558,7 @@ void PairAmoeba::multipole_real()
// resolve site torques then increment forces and virial
for (i = 0; i < nlocal; i++) {
torque2force(i,tq[i],fix,fiy,fiz,f);
torque2force(i,tq[i],fix,fiy,fiz,f);
iz = zaxis2local[i];
ix = xaxis2local[i];
@ -575,12 +575,12 @@ void PairAmoeba::multipole_real()
ziy = x[iy][2] - x[i][2];
vxx = xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
vxy = 0.5 * (yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
vxy = 0.5 * (yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
vxz = 0.5 * (zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
vxz = 0.5 * (zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
vyy = yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
vzz = zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
@ -643,10 +643,10 @@ void PairAmoeba::multipole_kspace()
double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
felec = electric / am_dielectric;
// perform dynamic allocation of arrays
// NOTE: just do this one time?
memory->create(cmp,nlocal,10,"ameoba/mpole:cmp");
memory->create(fmp,nlocal,10,"ameoba/mpole:fmp");
memory->create(cphi,nlocal,10,"ameoba/mpole:cphi");
@ -692,11 +692,11 @@ void PairAmoeba::multipole_kspace()
// pre-convolution operations including forward FFT
// gridfft = my portion of complex 3d grid in FFT decomp as 1d vector
double *gridfft = m_kspace->pre_convolution();
// ---------------------
// convolution operation
// convolution operation
// ---------------------
// zero virial accumulation variables
@ -718,7 +718,7 @@ void PairAmoeba::multipole_kspace()
pterm = pow((MY_PI/aewald),2.0);
volterm = MY_PI * volbox;
// printf("bsmod1 %e %e %e %e %e %e \n",bsmod1[0],bsmod1[1],bsmod1[2],bsmod1[3],bsmod1[4],bsmod1[5]);
// printf("bsmod2 %e %e %e %e %e %e \n",bsmod2[0],bsmod2[1],bsmod2[2],bsmod2[3],bsmod2[4],bsmod2[5]);
// printf("bsmod3 %e %e %e %e %e %e \n",bsmod3[0],bsmod3[1],bsmod3[2],bsmod3[3],bsmod3[4],bsmod3[5]);
@ -727,33 +727,33 @@ void PairAmoeba::multipole_kspace()
for (k = nzlo; k <= nzhi; k++) {
for (j = nylo; j <= nyhi; j++) {
for (i = nxlo; i <= nxhi; i++) {
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
r1 = (i >= nhalf1) ? i-nfft1 : i;
r2 = (j >= nhalf2) ? j-nfft2 : j;
r3 = (k >= nhalf3) ? k-nfft3 : k;
//printf("ijk %i %i %i r1r2r3 %f %f %f \n",i,j,k,r1,r2,r3);
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
term = -pterm * hsq;
expterm = 0.0;
if (term > -50.0 && hsq != 0.0) {
denom = volterm*hsq*bsmod1[i]*bsmod2[j]*bsmod3[k];
expterm = exp(term) / denom;
h1 = recip[0][0]*r1 + recip[0][1]*r2 + recip[0][2]*r3; // matvec
h2 = recip[1][0]*r1 + recip[1][1]*r2 + recip[1][2]*r3;
h3 = recip[2][0]*r1 + recip[2][1]*r2 + recip[2][2]*r3;
hsq = h1*h1 + h2*h2 + h3*h3;
term = -pterm * hsq;
expterm = 0.0;
if (term > -50.0 && hsq != 0.0) {
denom = volterm*hsq*bsmod1[i]*bsmod2[j]*bsmod3[k];
expterm = exp(term) / denom;
//printf("bsmod %e %e %e expterm %e \n",bsmod1[i],bsmod2[j],bsmod3[k],expterm);
struc2 = gridfft[n]*gridfft[n] + gridfft[n+1]*gridfft[n+1];
eterm = 0.5 * felec * expterm * struc2;
vterm = (2.0/hsq) * (1.0-term) * eterm;
vxx += h1*h1*vterm - eterm;
vyy += h2*h2*vterm - eterm;
vzz += h3*h3*vterm - eterm;
vxy += h1*h2*vterm;
vxz += h1*h3*vterm;
vyz += h2*h3*vterm;
}
struc2 = gridfft[n]*gridfft[n] + gridfft[n+1]*gridfft[n+1];
eterm = 0.5 * felec * expterm * struc2;
vterm = (2.0/hsq) * (1.0-term) * eterm;
vxx += h1*h1*vterm - eterm;
vyy += h2*h2*vterm - eterm;
vzz += h3*h3*vterm - eterm;
vxy += h1*h2*vterm;
vxz += h1*h3*vterm;
vyz += h2*h3*vterm;
}
gridfft[n] *= expterm;
gridfft[n+1] *= expterm;
n += 2;
n += 2;
}
}
}
@ -773,7 +773,7 @@ void PairAmoeba::multipole_kspace()
double ***gridpost = (double ***) m_kspace->post_convolution();
// get potential
fphi_mpole(gridpost,fphi);
//printf("fphi %e %e %e %e \n",fphi[0][0],fphi[0][1],fphi[0][2],fphi[0][3]);
@ -813,41 +813,41 @@ void PairAmoeba::multipole_kspace()
}
empole += 0.5*e;
//printf("mpole_force %g %g %g \n", f[0][0], f[0][1], f[0][2]);
// augment the permanent multipole virial contributions
for (i = 0; i < nlocal; i++) {
vxx = vxx - cmp[i][1]*cphi[i][1] - 2.0*cmp[i][4]*cphi[i][4] -
vxx = vxx - cmp[i][1]*cphi[i][1] - 2.0*cmp[i][4]*cphi[i][4] -
cmp[i][7]*cphi[i][7] - cmp[i][8]*cphi[i][8];
vxy = vxy - 0.5*(cmp[i][2]*cphi[i][1]+cmp[i][1]*cphi[i][2]) -
(cmp[i][4]+cmp[i][5])*cphi[i][7] - 0.5*cmp[i][7]*(cphi[i][4]+cphi[i][5]) -
vxy = vxy - 0.5*(cmp[i][2]*cphi[i][1]+cmp[i][1]*cphi[i][2]) -
(cmp[i][4]+cmp[i][5])*cphi[i][7] - 0.5*cmp[i][7]*(cphi[i][4]+cphi[i][5]) -
0.5*(cmp[i][8]*cphi[i][9]+cmp[i][9]*cphi[i][8]);
vxz = vxz - 0.5*(cmp[i][3]*cphi[i][1]+cmp[i][1]*cphi[i][3]) -
(cmp[i][4]+cmp[i][6])*cphi[i][8] - 0.5*cmp[i][8]*(cphi[i][4]+cphi[i][6]) -
vxz = vxz - 0.5*(cmp[i][3]*cphi[i][1]+cmp[i][1]*cphi[i][3]) -
(cmp[i][4]+cmp[i][6])*cphi[i][8] - 0.5*cmp[i][8]*(cphi[i][4]+cphi[i][6]) -
0.5*(cmp[i][7]*cphi[i][9]+cmp[i][9]*cphi[i][7]);
vyy = vyy - cmp[i][2]*cphi[i][2] - 2.0*cmp[i][5]*cphi[i][5] -
vyy = vyy - cmp[i][2]*cphi[i][2] - 2.0*cmp[i][5]*cphi[i][5] -
cmp[i][7]*cphi[i][7] - cmp[i][9]*cphi[i][9];
vyz = vyz - 0.5*(cmp[i][3]*cphi[i][2]+cmp[i][2]*cphi[i][3]) -
(cmp[i][5]+cmp[i][6])*cphi[i][9] - 0.5*cmp[i][9]*(cphi[i][5]+cphi[i][6]) -
vyz = vyz - 0.5*(cmp[i][3]*cphi[i][2]+cmp[i][2]*cphi[i][3]) -
(cmp[i][5]+cmp[i][6])*cphi[i][9] - 0.5*cmp[i][9]*(cphi[i][5]+cphi[i][6]) -
0.5*(cmp[i][7]*cphi[i][8]+cmp[i][8]*cphi[i][7]);
vzz = vzz - cmp[i][3]*cphi[i][3] - 2.0*cmp[i][6]*cphi[i][6] -
vzz = vzz - cmp[i][3]*cphi[i][3] - 2.0*cmp[i][6]*cphi[i][6] -
cmp[i][8]*cphi[i][8] - cmp[i][9]*cphi[i][9];
}
// resolve site torques then increment forces and virial
for (i = 0; i < nlocal; i++) {
tem[0] = cmp[i][3]*cphi[i][2] - cmp[i][2]*cphi[i][3] +
2.0*(cmp[i][6]-cmp[i][5])*cphi[i][9] +
cmp[i][8]*cphi[i][7] + cmp[i][9]*cphi[i][5] -
tem[0] = cmp[i][3]*cphi[i][2] - cmp[i][2]*cphi[i][3] +
2.0*(cmp[i][6]-cmp[i][5])*cphi[i][9] +
cmp[i][8]*cphi[i][7] + cmp[i][9]*cphi[i][5] -
cmp[i][7]*cphi[i][8] - cmp[i][9]*cphi[i][6];
tem[1] = cmp[i][1]*cphi[i][3] - cmp[i][3]*cphi[i][1] +
2.0*(cmp[i][4]-cmp[i][6])*cphi[i][8] +
cmp[i][7]*cphi[i][9] + cmp[i][8]*cphi[i][6] -
tem[1] = cmp[i][1]*cphi[i][3] - cmp[i][3]*cphi[i][1] +
2.0*(cmp[i][4]-cmp[i][6])*cphi[i][8] +
cmp[i][7]*cphi[i][9] + cmp[i][8]*cphi[i][6] -
cmp[i][8]*cphi[i][4] - cmp[i][9]*cphi[i][7];
tem[2] = cmp[i][2]*cphi[i][1] - cmp[i][1]*cphi[i][2] +
2.0*(cmp[i][5]-cmp[i][4])*cphi[i][7] +
cmp[i][7]*cphi[i][4] + cmp[i][9]*cphi[i][8] -
tem[2] = cmp[i][2]*cphi[i][1] - cmp[i][1]*cphi[i][2] +
2.0*(cmp[i][5]-cmp[i][4])*cphi[i][7] +
cmp[i][7]*cphi[i][4] + cmp[i][9]*cphi[i][8] -
cmp[i][7]*cphi[i][5] - cmp[i][8]*cphi[i][9];
torque2force(i,tem,fix,fiy,fiz,f);
@ -867,16 +867,16 @@ void PairAmoeba::multipole_kspace()
ziy = x[iy][2] - x[i][2];
vxx += xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
vxy += 0.5*(yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
vxy += 0.5*(yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
vxz += 0.5*(zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
vxz += 0.5*(zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
vyy += yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
vyz += 0.5*(zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
vyz += 0.5*(zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
vzz += zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
}
// increment total internal virial tensor components
virmpole[0] += vxx;
@ -947,7 +947,7 @@ void PairAmoeba::damppole(double r, int rorder, double alphai, double alphak,
dmpi[2] = 1.0 - (1.0 + dampi + 0.5*dampi2)*expi;
dmpi[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0)*expi;
dmpi[6] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 + dampi4/30.0)*expi;
dmpi[8] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dmpi[8] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
4.0*dampi4/105.0 + dampi5/210.0)*expi;
if (diff < eps) {
dmpk[0] = dmpi[0];
@ -964,7 +964,7 @@ void PairAmoeba::damppole(double r, int rorder, double alphai, double alphak,
dmpk[2] = 1.0 - (1.0 + dampk + 0.5*dampk2)*expk;
dmpk[4] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0)*expk;
dmpk[6] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + dampk4/30.0)*expk;
dmpk[8] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0 +
dmpk[8] = 1.0 - (1.0 + dampk + 0.5*dampk2 + dampk3/6.0 +
4.0*dampk4/105.0 + dampk5/210.0)*expk;
}
@ -973,21 +973,21 @@ void PairAmoeba::damppole(double r, int rorder, double alphai, double alphak,
if (diff < eps) {
dampi6 = dampi3 * dampi3;
dampi7 = dampi3 * dampi4;
dmpik[0] = 1.0 - (1.0 + 11.0*dampi/16.0 + 3.0*dampi2/16.0 +
dmpik[0] = 1.0 - (1.0 + 11.0*dampi/16.0 + 3.0*dampi2/16.0 +
dampi3/48.0)*expi;
dmpik[2] = 1.0 - (1.0 + dampi + 0.5*dampi2 +
dmpik[2] = 1.0 - (1.0 + dampi + 0.5*dampi2 +
7.0*dampi3/48.0 + dampi4/48.0)*expi;
dmpik[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dmpik[4] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/144.0)*expi;
dmpik[6] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dmpik[6] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/120.0 + dampi6/720.0)*expi;
dmpik[8] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/120.0 + dampi6/720.0 +
dmpik[8] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/120.0 + dampi6/720.0 +
dampi7/5040.0)*expi;
if (rorder >= 11) {
dampi8 = dampi4 * dampi4;
dmpik[10] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/120.0 + dampi6/720.0 +
dmpik[10] = 1.0 - (1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
dampi4/24.0 + dampi5/120.0 + dampi6/720.0 +
dampi7/5040.0 + dampi8/45360.0)*expi;
}
@ -998,41 +998,41 @@ void PairAmoeba::damppole(double r, int rorder, double alphai, double alphak,
termk = alphai2 / (alphai2-alphak2);
termi2 = termi * termi;
termk2 = termk * termk;
dmpik[0] = 1.0 - termi2*(1.0 + 2.0*termk + 0.5*dampi)*expi -
dmpik[0] = 1.0 - termi2*(1.0 + 2.0*termk + 0.5*dampi)*expi -
termk2*(1.0 + 2.0*termi + 0.5*dampk)*expk;
dmpik[2] = 1.0 - termi2*(1.0+dampi+0.5*dampi2)*expi -
termk2*(1.0+dampk+0.5*dampk2)*expk -
2.0*termi2*termk*(1.0+dampi)*expi -
2.0*termk2*termi*(1.0+dampk)*expk;
dmpik[4] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0)*expk -
2.0*termi2*termk*(1.0 + dampi + dampi2/3.0)*expi -
dmpik[4] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0)*expk -
2.0*termi2*termk*(1.0 + dampi + dampi2/3.0)*expi -
2.0*termk2*termi*(1.0 + dampk + dampk2/3.0)*expk;
dmpik[6] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 +
dampi3/6.0 + dampi4/30.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + dampk4/30.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 2.0*dampi2/5.0 + dampi3/15.0)*expi -
dmpik[6] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 +
dampi3/6.0 + dampi4/30.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + dampk4/30.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 2.0*dampi2/5.0 + dampi3/15.0)*expi -
2.0*termk2*termi*(1.0 + dampk + 2.0*dampk2/5.0 + dampk3/15.0)*expk;
dmpik[8] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
4.0*dampi4/105.0 + dampi5/210.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 +
4.0*dampk4/105.0 + dampk5/210.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 3.0*dampi2/7.0 +
2.0*dampi3/21.0 + dampi4/105.0)*expi -
2.0*termk2*termi*(1.0 + dampk + 3.0*dampk2/7.0 +
dmpik[8] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
4.0*dampi4/105.0 + dampi5/210.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 +
4.0*dampk4/105.0 + dampk5/210.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 3.0*dampi2/7.0 +
2.0*dampi3/21.0 + dampi4/105.0)*expi -
2.0*termk2*termi*(1.0 + dampk + 3.0*dampk2/7.0 +
2.0*dampk3/21.0 + dampk4/105.0)*expk;
if (rorder >= 11) {
dampi6 = dampi3 * dampi3;
dampk6 = dampk3 * dampk3;
dmpik[10] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
5.0*dampi4/126.0 + 2.0*dampi5/315.0 +
dampi6/1890.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + 5.0*dampk4/126.0 +
2.0*dampk5/315.0 + dampk6/1890.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 4.0*dampi2/9.0 + dampi3/9.0 +
dampi4/63.0 + dampi5/945.0)*expi -
2.0*termk2*termi*(1.0 + dampk + 4.0*dampk2/9.0 + dampk3/9.0 +
dmpik[10] = 1.0 - termi2*(1.0 + dampi + 0.5*dampi2 + dampi3/6.0 +
5.0*dampi4/126.0 + 2.0*dampi5/315.0 +
dampi6/1890.0)*expi -
termk2*(1.0 + dampk + 0.5*dampk2 + dampk3/6.0 + 5.0*dampk4/126.0 +
2.0*dampk5/315.0 + dampk6/1890.0)*expk -
2.0*termi2*termk*(1.0 + dampi + 4.0*dampi2/9.0 + dampi3/9.0 +
dampi4/63.0 + dampi5/945.0)*expi -
2.0*termk2*termi*(1.0 + dampk + 4.0*dampk2/9.0 + dampk3/9.0 +
dampk4/63.0 + dampk5/945.0)*expk;
}
}

500
src/AMOEBA/amoeba_polar.cpp
View File

File diff suppressed because it is too large Load Diff

114
src/AMOEBA/amoeba_repulsion.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -115,7 +115,7 @@ void PairAmoeba::repulsion()
// DEBUG
//FILE *fp = fopen("lammps.dat","w");
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = amtype[i];
@ -138,7 +138,7 @@ void PairAmoeba::repulsion()
qiyy = rpole[i][8];
qiyz = rpole[i][9];
qizz = rpole[i][12];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_repel = special_repel[sbmask15(j)];
@ -231,11 +231,11 @@ void PairAmoeba::repulsion()
dkqirx = dkqiz*yr - dkqiy*zr;
dkqiry = dkqix*zr - dkqiz*xr;
dkqirz = dkqiy*xr - dkqix*yr;
dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy -
dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy -
2.0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz-qixz*qkxy-qiyz*qkyy-qizz*qkyz);
dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz -
dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz -
2.0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz-qixx*qkxz-qixy*qkyz-qixz*qkzz);
dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix -
dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix -
2.0*(qixx*qkxy+qixy*qkyy+qixz*qkyz-qixy*qkxx-qiyy*qkxy-qiyz*qkxz);
// get reciprocal distance terms for this interaction
@ -258,7 +258,7 @@ void PairAmoeba::repulsion()
term3 = vali*qkr + valk*qir - dir*dkr + 2.0*(dkqi-diqk+qiqk);
term4 = dir*qkr - dkr*qir - 4.0*qik;
term5 = qir*qkr;
eterm = term1*dmpik[0] + term2*dmpik[2] +
eterm = term1*dmpik[0] + term2*dmpik[2] +
term3*dmpik[4] + term4*dmpik[6] + term5*dmpik[8];
// compute the Pauli repulsion energy for this interaction
@ -268,7 +268,7 @@ void PairAmoeba::repulsion()
// calculate intermediate terms for force and torque
de = term1*dmpik[2] + term2*dmpik[4] + term3*dmpik[6] +
de = term1*dmpik[2] + term2*dmpik[4] + term3*dmpik[6] +
term4*dmpik[8] + term5*dmpik[10];
term1 = -valk*dmpik[2] + dkr*dmpik[4] - qkr*dmpik[6];
term2 = vali*dmpik[2] + dir*dmpik[4] + qir*dmpik[6];
@ -279,11 +279,11 @@ void PairAmoeba::repulsion()
// compute the force components for this interaction
frcx = de*xr + term1*dix + term2*dkx + term3*(diqkx-dkqix) +
frcx = de*xr + term1*dix + term2*dkx + term3*(diqkx-dkqix) +
term4*qix + term5*qkx + term6*(qixk+qkxi);
frcy = de*yr + term1*diy + term2*dky + term3*(diqky-dkqiy) +
frcy = de*yr + term1*diy + term2*dky + term3*(diqky-dkqiy) +
term4*qiy + term5*qky + term6*(qiyk+qkyi);
frcz = de*zr + term1*diz + term2*dkz + term3*(diqkz-dkqiz) +
frcz = de*zr + term1*diz + term2*dkz + term3*(diqkz-dkqiz) +
term4*qiz + term5*qkz + term6*(qizk+qkzi);
frcx = frcx*rr1 + eterm*rr3*xr;
frcy = frcy*rr1 + eterm*rr3*yr;
@ -293,18 +293,18 @@ void PairAmoeba::repulsion()
frcz = sizik * frcz;
// compute the torque components for this interaction
ttri[0] = -dmpik[2]*dikx + term1*dirx + term3*(dqikx+dkqirx) -
ttri[0] = -dmpik[2]*dikx + term1*dirx + term3*(dqikx+dkqirx) -
term4*qirx - term6*(qikrx+qikx);
ttri[1] = -dmpik[2]*diky + term1*diry + term3*(dqiky+dkqiry) -
ttri[1] = -dmpik[2]*diky + term1*diry + term3*(dqiky+dkqiry) -
term4*qiry - term6*(qikry+qiky);
ttri[2] = -dmpik[2]*dikz + term1*dirz + term3*(dqikz+dkqirz) -
ttri[2] = -dmpik[2]*dikz + term1*dirz + term3*(dqikz+dkqirz) -
term4*qirz - term6*(qikrz+qikz);
ttrk[0] = dmpik[2]*dikx + term2*dkrx - term3*(dqikx+diqkrx) -
ttrk[0] = dmpik[2]*dikx + term2*dkrx - term3*(dqikx+diqkrx) -
term5*qkrx - term6*(qkirx-qikx);
ttrk[1] = dmpik[2]*diky + term2*dkry - term3*(dqiky+diqkry) -
ttrk[1] = dmpik[2]*diky + term2*dkry - term3*(dqiky+diqkry) -
term5*qkry - term6*(qkiry-qiky);
ttrk[2] = dmpik[2]*dikz + term2*dkrz - term3*(dqikz+diqkrz) -
ttrk[2] = dmpik[2]*dikz + term2*dkrz - term3*(dqikz+diqkrz) -
term5*qkrz - term6*(qkirz-qikz);
ttri[0] = sizik * ttri[0] * rr1;
ttri[1] = sizik * ttri[1] * rr1;
@ -371,7 +371,7 @@ void PairAmoeba::repulsion()
// energy = e
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}
@ -379,7 +379,7 @@ void PairAmoeba::repulsion()
// DEBUG
//fclose(fp);
// reverse comm to sum torque from ghost atoms to owned atoms
crstyle = TORQUE;
@ -407,11 +407,11 @@ void PairAmoeba::repulsion()
vxx = xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
vyy = yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
vzz = zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
vxy = 0.5 * (yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
vxy = 0.5 * (yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
vxz = 0.5 * (zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
vxz = 0.5 * (zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
virrepulse[0] += vxx;
@ -423,7 +423,7 @@ void PairAmoeba::repulsion()
// virial = 6-vec vir
// NOTE: add tally function
if (evflag) {
}
}
@ -440,7 +440,7 @@ void PairAmoeba::repulsion()
150, 084104 (2019)
------------------------------------------------------------------------- */
void PairAmoeba::damprep(double r, double r2, double rr1, double rr3,
void PairAmoeba::damprep(double r, double r2, double rr1, double rr3,
double rr5, double rr7, double rr9, double rr11,
int rorder, double dmpi, double dmpk, double *dmpik)
{
@ -519,46 +519,46 @@ void PairAmoeba::damprep(double r, double r2, double rr1, double rr3,
tmp = 4.0 * dmpi2 * dmpk2 / term;
s = (dampi-tmp)*expk + (dampk+tmp)*expi;
ds = (dmpi2*dmpk2*r2 - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
ds = (dmpi2*dmpk2*r2 - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi2*dmpk2*r2 + 4.0*dmpi22*dmpk2*r/term + 4.0*dmpi2*dmpk2/term) * expi;
d2s = (dmpi2*dmpk2*r2/3.0 + dmpi2*dmpk22*r3/3.0 -
(4.0/3.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi2*dmpk2*r2/3.0 + dmpi22*dmpk2*r3/3.0 +
(4.0/3.0)*dmpi23*dmpk2*r2/term + 4.0*dmpi22*dmpk2*r/term +
d2s = (dmpi2*dmpk2*r2/3.0 + dmpi2*dmpk22*r3/3.0 -
(4.0/3.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi2*dmpk2*r2/3.0 + dmpi22*dmpk2*r3/3.0 +
(4.0/3.0)*dmpi23*dmpk2*r2/term + 4.0*dmpi22*dmpk2*r/term +
4.0*dmpi2*dmpk2/term) * expi;
d3s = (dmpi2*dmpk23*r4/15.0 + dmpi2*dmpk22*r3/5.0 + dmpi2*dmpk2*r2/5.0 -
(4.0/15.0)*dmpi2*dmpk24*r3/term - (8.0/5.0)*dmpi2*dmpk23*r2/term -
4.0*dmpi2*dmpk22*r/term - 4.0/term*dmpi2*dmpk2) * expk +
(dmpi23*dmpk2*r4/15.0 + dmpi22*dmpk2*r3/5.0 + dmpi2*dmpk2*r2/5.0 +
(4.0/15.0)*dmpi24*dmpk2*r3/term + (8.0/5.0)*dmpi23*dmpk2*r2/term +
d3s = (dmpi2*dmpk23*r4/15.0 + dmpi2*dmpk22*r3/5.0 + dmpi2*dmpk2*r2/5.0 -
(4.0/15.0)*dmpi2*dmpk24*r3/term - (8.0/5.0)*dmpi2*dmpk23*r2/term -
4.0*dmpi2*dmpk22*r/term - 4.0/term*dmpi2*dmpk2) * expk +
(dmpi23*dmpk2*r4/15.0 + dmpi22*dmpk2*r3/5.0 + dmpi2*dmpk2*r2/5.0 +
(4.0/15.0)*dmpi24*dmpk2*r3/term + (8.0/5.0)*dmpi23*dmpk2*r2/term +
4.0*dmpi22*dmpk2*r/term + 4.0/term*dmpi2*dmpk2) * expi;
d4s = (dmpi2*dmpk24*r5/105.0 + (2.0/35.0)*dmpi2*dmpk23*r4 +
dmpi2*dmpk22*r3/7.0 + dmpi2*dmpk2*r2/7.0 -
(4.0/105.0)*dmpi2*dmpk25*r4/term - (8.0/21.0)*dmpi2*dmpk24*r3/term -
(12.0/7.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi24*dmpk2*r5/105.0 + (2.0/35.0)*dmpi23*dmpk2*r4 +
dmpi22*dmpk2*r3/7.0 + dmpi2*dmpk2*r2/7.0 +
(4.0/105.0)*dmpi25*dmpk2*r4/term + (8.0/21.0)*dmpi24*dmpk2*r3/term +
(12.0/7.0)*dmpi23*dmpk2*r2/term + 4.0*dmpi22*dmpk2*r/term +
d4s = (dmpi2*dmpk24*r5/105.0 + (2.0/35.0)*dmpi2*dmpk23*r4 +
dmpi2*dmpk22*r3/7.0 + dmpi2*dmpk2*r2/7.0 -
(4.0/105.0)*dmpi2*dmpk25*r4/term - (8.0/21.0)*dmpi2*dmpk24*r3/term -
(12.0/7.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi24*dmpk2*r5/105.0 + (2.0/35.0)*dmpi23*dmpk2*r4 +
dmpi22*dmpk2*r3/7.0 + dmpi2*dmpk2*r2/7.0 +
(4.0/105.0)*dmpi25*dmpk2*r4/term + (8.0/21.0)*dmpi24*dmpk2*r3/term +
(12.0/7.0)*dmpi23*dmpk2*r2/term + 4.0*dmpi22*dmpk2*r/term +
4.0*dmpi2*dmpk2/term) * expi;
if (rorder >= 11) {
r6 = r5 * r;
dmpi26 = dmpi25 * dmpi2;
dmpk26 = dmpk25 * dmpk2;
d5s = (dmpi2*dmpk25*r6/945.0 + (2.0/189.0)*dmpi2*dmpk24*r5 +
dmpi2*dmpk23*r4/21.0 + dmpi2*dmpk22*r3/9.0 + dmpi2*dmpk2*r2/9.0 -
(4.0/945.0)*dmpi2*dmpk26*r5/term -
(4.0/63.0)*dmpi2*dmpk25*r4/term - (4.0/9.0)*dmpi2*dmpk24*r3/term -
(16.0/9.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi25*dmpk2*r6/945.0 + (2.0/189.0)*dmpi24*dmpk2*r5 +
dmpi23*dmpk2*r4/21.0 + dmpi22*dmpk2*r3/9.0 + dmpi2*dmpk2*r2/9.0 +
(4.0/945.0)*dmpi26*dmpk2*r5/term + (4.0/63.0)*dmpi25*dmpk2*r4/term +
(4.0/9.0)*dmpi24*dmpk2*r3/term + (16.0/9.0)*dmpi23*dmpk2*r2/term +
d5s = (dmpi2*dmpk25*r6/945.0 + (2.0/189.0)*dmpi2*dmpk24*r5 +
dmpi2*dmpk23*r4/21.0 + dmpi2*dmpk22*r3/9.0 + dmpi2*dmpk2*r2/9.0 -
(4.0/945.0)*dmpi2*dmpk26*r5/term -
(4.0/63.0)*dmpi2*dmpk25*r4/term - (4.0/9.0)*dmpi2*dmpk24*r3/term -
(16.0/9.0)*dmpi2*dmpk23*r2/term - 4.0*dmpi2*dmpk22*r/term -
4.0*dmpi2*dmpk2/term) * expk +
(dmpi25*dmpk2*r6/945.0 + (2.0/189.0)*dmpi24*dmpk2*r5 +
dmpi23*dmpk2*r4/21.0 + dmpi22*dmpk2*r3/9.0 + dmpi2*dmpk2*r2/9.0 +
(4.0/945.0)*dmpi26*dmpk2*r5/term + (4.0/63.0)*dmpi25*dmpk2*r4/term +
(4.0/9.0)*dmpi24*dmpk2*r3/term + (16.0/9.0)*dmpi23*dmpk2*r2/term +
4.0*dmpi22*dmpk2*r/term + 4.0*dmpi2*dmpk2/term) * expi;
}
}

230
src/AMOEBA/amoeba_utils.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -59,7 +59,7 @@ void PairAmoeba::kmpole()
double **pole = fixpole->astore;
int **nspecial = atom->nspecial;
int **special = atom->special;
tagint **special = atom->special;
int nlocal = atom->nlocal;
int nmissing = 0;
@ -111,49 +111,49 @@ void PairAmoeba::kmpole()
ztype = zpole[itype][iframe];
for (j12 = 0; j12 < nspecial[i][0]; j12++) {
jneigh = bondneigh[j12];
j = atom->map(jneigh);
if (j < 0)
j = atom->map(jneigh);
if (j < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
jtype = amtype[j];
if (jtype == ztype) {
for (k12 = 0; k12 < nspecial[i][0]; k12++) {
if (k12 == j12) continue;
jtype = amtype[j];
if (jtype == ztype) {
for (k12 = 0; k12 < nspecial[i][0]; k12++) {
if (k12 == j12) continue;
kneigh = bondneigh[k12];
k = atom->map(kneigh);
if (k < 0)
k = atom->map(kneigh);
if (k < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
ktype = amtype[k];
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
} else {
for (m12 = 0; m12 < nspecial[i][0]; m12++) {
if (m12 == j12 || m12 == k12) continue;
ktype = amtype[k];
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
} else {
for (m12 = 0; m12 < nspecial[i][0]; m12++) {
if (m12 == j12 || m12 == k12) continue;
mneigh = bondneigh[m12];
m = atom->map(mneigh);
if (m < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
mtype = amtype[m];
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
}
}
}
}
m = atom->map(mneigh);
if (m < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
mtype = amtype[m];
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
}
}
}
}
}
}
@ -167,55 +167,55 @@ void PairAmoeba::kmpole()
ztype = zpole[itype][iframe];
for (j12 = 0; j12 < nspecial[i][0]; j12++) {
jneigh = bondneigh[j12];
j = atom->map(jneigh);
if (j < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
jtype = amtype[j];
if (jtype == ztype) {
for (k13 = nspecial[i][0]; k13 < nspecial[i][1]; k13++) {
j = atom->map(jneigh);
if (j < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
jtype = amtype[j];
if (jtype == ztype) {
for (k13 = nspecial[i][0]; k13 < nspecial[i][1]; k13++) {
kneigh = angleneigh[k13];
k = atom->map(kneigh);
if (k < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
ktype = amtype[k];
path = false;
for (m12 = 0; m12 < nspecial[k][0]; m12++)
if (special[k][m12] == jneigh) path = true;
if (!path) continue;
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
} else {
for (m13 = nspecial[i][0]; m13 < nspecial[i][1]; m13++) {
if (m13 == k13) continue;
k = atom->map(kneigh);
if (k < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
ktype = amtype[k];
path = false;
for (m12 = 0; m12 < nspecial[k][0]; m12++)
if (special[k][m12] == jneigh) path = true;
if (!path) continue;
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
} else {
for (m13 = nspecial[i][0]; m13 < nspecial[i][1]; m13++) {
if (m13 == k13) continue;
mneigh = angleneigh[m13];
m = atom->map(mneigh);
if (m < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
mtype = amtype[m];
path = false;
for (m12 = 0; m12 < nspecial[m][0]; m12++)
if (special[m][m12] == jneigh) path = true;
if (!path) continue;
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
}
}
}
}
m = atom->map(mneigh);
if (m < 0)
error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
mtype = amtype[m];
path = false;
for (m12 = 0; m12 < nspecial[m][0]; m12++)
if (special[m][m12] == jneigh) path = true;
if (!path) continue;
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
}
}
}
}
}
}
@ -229,19 +229,19 @@ void PairAmoeba::kmpole()
ztype = zpole[itype][iframe];
for (j12 = 0; j12 < nspecial[i][0]; j12++) {
jneigh = bondneigh[j12];
j = atom->map(jneigh);
if (j < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
jtype = amtype[j];
if (jtype == ztype) {
if (xtype == 0 && !flag) {
flag = 1;
zaxis[i] = jtype;
xaxis[i] = yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
j = atom->map(jneigh);
if (j < 0) error->one(FLERR,"AMOEBA kmpole() could not find bond partner");
jtype = amtype[j];
if (jtype == ztype) {
if (xtype == 0 && !flag) {
flag = 1;
zaxis[i] = jtype;
xaxis[i] = yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
}
}
@ -254,11 +254,11 @@ void PairAmoeba::kmpole()
ytype = ypole[itype][iframe];
ztype = zpole[itype][iframe];
if (ztype == 0 && !flag) {
flag = 1;
zaxis[i] = xaxis[i] = yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
flag = 1;
zaxis[i] = xaxis[i] = yaxis[i] = 0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
}
}
@ -639,7 +639,7 @@ void PairAmoeba::add_onefive_neighbors()
tagint *tag = atom->tag;
int *nspecial15 = atom->nspecial15;
int **special15 = atom->special15;
tagint **special15 = atom->special15;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@ -647,7 +647,7 @@ void PairAmoeba::add_onefive_neighbors()
list15 = special15[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
which = j >> SBBITS & 3;
@ -655,12 +655,12 @@ void PairAmoeba::add_onefive_neighbors()
jtag = tag[j];
if (!which) {
for (k = 0; k < n15; k++) {
if (list15[k] == jtag) {
which = 4;
break;
}
}
for (k = 0; k < n15; k++) {
if (list15[k] == jtag) {
which = 4;
break;
}
}
}
if (which) jlist[jj] = j ^ (which << SBBITS15);
@ -792,7 +792,7 @@ void PairAmoeba::find_multipole_neighbors()
f = force vector for local atoms, multiple atoms will be updated
------------------------------------------------------------------------- */
void PairAmoeba::torque2force(int i, double *trq,
void PairAmoeba::torque2force(int i, double *trq,
double *frcx, double *frcy, double *frcz,
double **f)
{

42
src/AMOEBA/angle_amoeba.cpp
View File

@ -33,7 +33,7 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
{
pflag = nullptr;
ubflag = nullptr;
@ -122,11 +122,11 @@ void AngleAmoeba::compute(int eflag, int vflag)
tinker_angle(i1,i2,i3,type,eflag);
// bondangle = bond-stretch cross term in Tinker
if (ba_k1[type] != 0.0)
tinker_bondangle(i1,i2,i3,type,eflag);
}
// Urey-Bradley H-H bond term within water molecules
if (enable_urey) {
@ -162,7 +162,7 @@ void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
@ -171,7 +171,7 @@ void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
@ -181,7 +181,7 @@ void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL;
s = 1.0/s;
// force & energy for angle term
dtheta = acos(c) - theta0[type];
@ -198,17 +198,17 @@ void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
a22 = a*c / rsq2;
f1[0] = a11*delx1 + a12*delx2;
f1[1] = a11*dely1 + a12*dely2;
f1[2] = a11*delz1 + a12*delz2;
f3[0] = a22*delx2 + a12*delx1;
f3[1] = a22*dely2 + a12*dely1;
f3[2] = a22*delz2 + a12*delz1;
eangle = 0.0;
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
// apply force to each of 3 atoms
@ -261,7 +261,7 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
int newton_bond = force->newton_bond;
// i4 = index of third atom that i2 is bonded to
i1tag = atom->tag[i1];
i3tag = atom->tag[i3];
@ -329,7 +329,7 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
cosine = MIN(1.0,MAX(-1.0,cosine));
// force & energy for angle term
dtheta = acos(cosine) - theta0[type];
dtheta2 = dtheta*dtheta;
dtheta3 = dtheta2*dtheta;
@ -341,7 +341,7 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
4.0*k4[type]*dtheta3 + 5.0*k5[type]*dtheta4 + 6.0*k6[type]*dtheta5;
eangle = 0.0;
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 +
k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
// chain rule terms for first derivative components
@ -445,7 +445,7 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
@ -454,7 +454,7 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
@ -503,7 +503,7 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
f1[0] = -(vx11 + b1*delx1 + vx12);
f1[1] = -(vy11 + b1*dely1 + vy12);
f1[2] = -(vz11 + b1*delz1 + vz12);
f3[0] = -(vx21 + b2*delx2 + vx22);
f3[1] = -(vy21 + b2*dely2 + vy22);
f3[2] = -(vz21 + b2*delz2 + vz22);
@ -558,7 +558,7 @@ void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
rk = ub_k[type] * dr;
// force & energy
if (r > 0.0) fbond = -2.0*rk/r;
else fbond = 0.0;
@ -610,7 +610,7 @@ void AngleAmoeba::allocate()
memory->create(setflag_ba,n+1,"angle:setflag_ba");
memory->create(setflag_ub,n+1,"angle:setflag_ub");
for (int i = 1; i <= n; i++)
for (int i = 1; i <= n; i++)
setflag[i] = setflag_a[i] = setflag_ba[i] = setflag_ub[i] = 0;
}
@ -627,7 +627,7 @@ void AngleAmoeba::coeff(int narg, char **arg)
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
int count = 0;
if (strcmp(arg[1],"ba") == 0) {
if (narg != 6) error->all(FLERR,"Incorrect args for angle coefficients");
@ -658,7 +658,7 @@ void AngleAmoeba::coeff(int narg, char **arg)
count++;
}
} else {
} else {
if (narg != 9) error->all(FLERR,"Incorrect args for angle coefficients");
int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
@ -671,7 +671,7 @@ void AngleAmoeba::coeff(int narg, char **arg)
double k6_one = utils::numeric(FLERR,arg[8],false,lmp);
// convert theta0 from degrees to radians
for (int i = ilo; i <= ihi; i++) {
pflag[i] = pflag_one;
ubflag[i] = ubflag_one;
@ -689,7 +689,7 @@ void AngleAmoeba::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
for (int i = ilo; i <= ihi; i++)
if (setflag_a[i] == 1 && setflag_ba[i] == 1 && setflag_ub[i])
if (setflag_a[i] == 1 && setflag_ba[i] == 1 && setflag_ub[i])
setflag[i] = 1;
}

4
src/AMOEBA/atom_vec_amoeba.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -26,7 +26,7 @@ AtomVecAmoeba::AtomVecAmoeba(LAMMPS *lmp) : AtomVec(lmp)
atom->molecule_flag = atom->q_flag = 1;
atom->nspecial15_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter

2
src/AMOEBA/atom_vec_amoeba.h
View File

@ -1,6 +1,6 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract

76
src/AMOEBA/fix_amoeba_bitorsion.cpp
View File

@ -256,7 +256,7 @@ void FixAmoebaBiTorsion::pre_neighbor()
if (max_bitorsion_list == 0) {
if (nprocs == 1) max_bitorsion_list = nbitorsions;
else max_bitorsion_list =
else max_bitorsion_list =
static_cast<int> (LB_FACTOR*nbitorsions/nprocs);
memory->create(bitorsion_list,max_bitorsion_list,6,
"bitorsion:bitorsion_list");
@ -559,7 +559,7 @@ void FixAmoebaBiTorsion::post_force(int vflag)
dedxu = -dedang1 * (yu*zcb - ycb*zu) / (ru2*rcb);
dedyu = -dedang1 * (zu*xcb - zcb*xu) / (ru2*rcb);
dedzu = -dedang1 * (xu*ycb - xcb*yu) / (ru2*rcb);
// compute first derivative components for first angle
dedxia = zcb*dedyt - ycb*dedzt;
@ -576,18 +576,18 @@ void FixAmoebaBiTorsion::post_force(int vflag)
dedzid = ycb*dedxu - xcb*dedyu;
// chain rule terms for second angle derivative components
xec = xie - xic;
yec = yie - yic;
zec = zie - zic;
dedxu2 = dedang2 * (yu*zdc - ydc*zu) / (ru2*rdc);
dedyu2 = dedang2 * (zu*xdc - zdc*xu) / (ru2*rdc);
dedzu2 = dedang2 * (xu*ydc - xdc*yu) / (ru2*rdc);
dedxv2 = -dedang2 * (yv*zdc - ydc*zv) / (rv2*rdc);
dedyv2 = -dedang2 * (zv*xdc - zdc*xv) / (rv2*rdc);
dedzv2 = -dedang2 * (xv*ydc - xdc*yv) / (rv2*rdc);
// compute first derivative components for second angle
dedxib2 = zdc*dedyu2 - ydc*dedzu2;
@ -662,14 +662,14 @@ void FixAmoebaBiTorsion::post_force(int vflag)
vyy2 = ydc*(dedyid2+dedyie2) - ycb*dedyib2 + yed*dedyie2;
vzy2 = zdc*(dedyid2+dedyie2) - zcb*dedyib2 + zed*dedyie2;
vzz2 = zdc*(dedzid2+dedzie2) - zcb*dedzib2 + zed*dedzie2;
v[0] += vxx + vxx2;
v[1] += vyy + vyy2;
v[2] += vzz + vzz2;
v[3] += vyx + vyx2;
v[4] += vzx + vzx2;
v[5] += vzy + vzy2;
ev_tally(nlist,list,5.0,e,v);
}
}
@ -735,7 +735,7 @@ void FixAmoebaBiTorsion::read_grid_data(char *bitorsion_file)
eof = fgets(line,MAXLINE,fp);
eof = fgets(line,MAXLINE,fp);
eof = fgets(line,MAXLINE,fp);
if (eof == nullptr)
if (eof == nullptr)
error->one(FLERR,"Unexpected end of fix amoeba/bitorsion file");
sscanf(line,"%d",&nbitypes);
@ -762,7 +762,7 @@ void FixAmoebaBiTorsion::read_grid_data(char *bitorsion_file)
if (me == 0) {
eof = fgets(line,MAXLINE,fp);
eof = fgets(line,MAXLINE,fp);
if (eof == nullptr)
if (eof == nullptr)
error->one(FLERR,"Unexpected end of fix amoeba/bitorsion file");
sscanf(line,"%d %d %d",&tmp,&nx,&ny);
}
@ -782,7 +782,7 @@ void FixAmoebaBiTorsion::read_grid_data(char *bitorsion_file)
for (int iy = 0; iy < ny; iy++) {
for (int ix = 0; ix < nx; ix++) {
eof = fgets(line,MAXLINE,fp);
if (eof == nullptr)
if (eof == nullptr)
error->one(FLERR,"Unexpected end of fix amoeba/bitorsion file");
sscanf(line,"%lg %lg %lg",&xgrid,&ygrid,&value);
if (iy == 0) ttx[itype][ix] = xgrid;
@ -857,7 +857,7 @@ void FixAmoebaBiTorsion::create_splines()
// cyclic = 1 if angle range is -180.0 to 180.0 in both dims
// error if cyclic and x,y pairs of edges of 2d array values do not match
// equality comparisons are within eps
int cyclic = 1;
double eps = 1.0e-6;
@ -876,19 +876,19 @@ void FixAmoebaBiTorsion::create_splines()
// spline fit of derivatives about 1st torsion
for (int i = 0; i < nx; i++)
for (int i = 0; i < nx; i++)
tmp1[i] = ttx[itype][i];
for (int j = 0; j < ny; j++) {
for (int i = 0; i < nx; i++)
for (int i = 0; i < nx; i++)
tmp2[i] = tbf[itype][j*nx+i];
if (cyclic)
if (cyclic)
cspline(nx-1,tmp1,tmp2,bs,cs,ds,tmp3,tmp4,tmp5,tmp6,tmp7);
else
else
nspline(nx-1,tmp1,tmp2,bs,cs,tmp3,tmp4,tmp5,tmp6,tmp7);
for (int i = 0; i < nx; i++)
for (int i = 0; i < nx; i++)
tbx[itype][j*nx+i] = bs[i];
}
@ -901,9 +901,9 @@ void FixAmoebaBiTorsion::create_splines()
for (int j = 0; j < ny; j++)
tmp2[j] = tbf[itype][j*nx+i];
if (cyclic)
if (cyclic)
cspline(ny-1,tmp1,tmp2,bs,cs,ds,tmp3,tmp4,tmp5,tmp6,tmp7);
else
else
nspline(ny-1,tmp1,tmp2,bs,cs,tmp3,tmp4,tmp5,tmp6,tmp7);
for (int j = 0; j < ny; j++)
@ -919,12 +919,12 @@ void FixAmoebaBiTorsion::create_splines()
for (int j = 0; j < ny; j++)
tmp2[j] = tbx[itype][j*nx+i];
if (cyclic)
if (cyclic)
cspline(ny-1,tmp1,tmp2,bs,cs,ds,tmp3,tmp4,tmp5,tmp6,tmp7);
else
else
nspline(ny-1,tmp1,tmp2,bs,cs,tmp3,tmp4,tmp5,tmp6,tmp7);
for (int j = 0; j < ny; j++)
for (int j = 0; j < ny; j++)
tbxy[itype][j*nx+i] = bs[j];
}
}
@ -953,9 +953,9 @@ void FixAmoebaBiTorsion::create_splines()
rest of args are outputs
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::nspline(int n, double *x0, double *y0,
void FixAmoebaBiTorsion::nspline(int n, double *x0, double *y0,
double *s1, double *s2,
double *h, double *g, double *dy,
double *h, double *g, double *dy,
double *dla, double *dmu)
{
int i;
@ -982,7 +982,7 @@ void FixAmoebaBiTorsion::nspline(int n, double *x0, double *y0,
}
// set the initial value via natural boundary condition
dla[n] = 1.0;
dla[0] = 0.0;
dmu[n] = 0.0;
@ -1024,9 +1024,9 @@ void FixAmoebaBiTorsion::nspline(int n, double *x0, double *y0,
rest of args are outputs
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::cspline(int n, double *xn, double *fn,
void FixAmoebaBiTorsion::cspline(int n, double *xn, double *fn,
double *b, double *c, double *d,
double *h, double *du, double *dm,
double *h, double *du, double *dm,
double *rc, double *rs)
{
int i;
@ -1038,7 +1038,7 @@ void FixAmoebaBiTorsion::cspline(int n, double *xn, double *fn,
// get auxiliary variables and matrix elements on first call
for (i = 0; i < n; i++)
for (i = 0; i < n; i++)
h[i] = xn[i+1] - xn[i];
h[n] = h[0];
@ -1083,7 +1083,7 @@ void FixAmoebaBiTorsion::cspline(int n, double *xn, double *fn,
du,cr,rs,x,iflag are outputs
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::cytsy(int n, double *dm, double *du,
void FixAmoebaBiTorsion::cytsy(int n, double *dm, double *du,
double *cr, double *rs, double *x, int &iflag)
{
iflag = -2;
@ -1103,14 +1103,14 @@ void FixAmoebaBiTorsion::cytsy(int n, double *dm, double *du,
du,cr,iflag are outputs
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::cytsyp(int n, double *dm, double *du,
void FixAmoebaBiTorsion::cytsyp(int n, double *dm, double *du,
double *cr, int &iflag)
{
int i;
double temp1,temp2;
// set error bound and test for condition n greater than 2
double eps = 1.0e-8;
iflag = -2;
@ -1135,7 +1135,7 @@ void FixAmoebaBiTorsion::cytsyp(int n, double *dm, double *du,
return;
}
// factoring a while checking for a positive definite and
// factoring a while checking for a positive definite and
// strong nonsingular matrix a
temp1 = du[1];
@ -1177,7 +1177,7 @@ void FixAmoebaBiTorsion::cytsyp(int n, double *dm, double *du,
dm[n] = dm[n] - dm[n-1]*du[n-1]*du[n-1];
temp1 = 0.0;
for (i = 1; i < n-1; i++)
for (i = 1; i < n-1; i++)
temp1 += dm[i]*cr[i]*cr[i];
dm[n] = dm[n] - temp1;
@ -1201,7 +1201,7 @@ void FixAmoebaBiTorsion::cytsyp(int n, double *dm, double *du,
rs,x are outputs
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::cytsys(int n, double *dm, double *du,
void FixAmoebaBiTorsion::cytsys(int n, double *dm, double *du,
double *cr, double *rs, double *x)
{
int i;
@ -1233,7 +1233,7 @@ void FixAmoebaBiTorsion::cytsys(int n, double *dm, double *du,
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::chkttor(int ib, int ic, int id,
void FixAmoebaBiTorsion::chkttor(int ib, int ic, int id,
double &sign, double &value1, double &value2)
{
}
@ -1245,10 +1245,10 @@ void FixAmoebaBiTorsion::chkttor(int ib, int ic, int id,
output = ansy,ansy1,ansy2
------------------------------------------------------------------------- */
void FixAmoebaBiTorsion::bcuint1(double *y, double *y1,
void FixAmoebaBiTorsion::bcuint1(double *y, double *y1,
double *y2, double *y12,
double x1l, double x1u, double x2l, double x2u,
double x1, double x2,
double x1l, double x1u, double x2l, double x2u,
double x1, double x2,
double &ansy, double &ansy1, double &ansy2)
{
double c[4][4];
@ -1466,7 +1466,7 @@ void FixAmoebaBiTorsion::write_data_section_size(int /*mth*/, int &nx, int &ny)
for (i = 0; i < nlocal; i++)
for (m = 0; m < num_bitorsion[i]; m++)
if (bitorsion_atom3[i][m] == tag[i]) nx++;
ny = 6;
}

4
src/AMOEBA/fix_amoeba_bitorsion.h
View File

@ -111,9 +111,9 @@ class FixAmoebaBiTorsion : public Fix {
void chkttor(int, int, int, double &, double &, double &);
void bcuint1(double *, double *, double *, double *,
double, double, double, double, double, double,
double, double, double, double, double, double,
double &, double &, double &);
void bcucof(double *, double *, double *, double *, double, double,
void bcucof(double *, double *, double *, double *, double, double,
double [][4]);
};

24
src/AMOEBA/fix_amoeba_pitorsion.cpp
View File

@ -222,7 +222,7 @@ void FixAmoebaPiTorsion::pre_neighbor()
if (max_pitorsion_list == 0) {
if (nprocs == 1) max_pitorsion_list = npitorsions;
else max_pitorsion_list =
else max_pitorsion_list =
static_cast<int> (LB_FACTOR*npitorsions/nprocs);
memory->create(pitorsion_list,max_pitorsion_list,7,
"pitorsion:pitorsion_list");
@ -423,7 +423,7 @@ void FixAmoebaPiTorsion::post_force(int vflag)
sine = (xdc*xtu + ydc*ytu + zdc*ztu) / (rdc*rtru);
// set the pi-system torsion parameters for this angle
v2 = kpit[ptype];
c2 = -1.0;
s2 = 0.0;
@ -625,11 +625,11 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
sscanf(keyword,"%d",&npitorsion_types);
which = 1;
} else error->all(FLERR,"Invalid read data section for fix amoeba/pitorsion");
// loop over lines of PiTorsion Coeffs
// tokenize the line into values
// initialize kpit vector
if (which == 1) {
int itype;
double value;
@ -662,7 +662,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
*next = '\0';
int nwords = utils::count_words(utils::trim_comment(buf));
*next = '\n';
if (nwords != 8)
error->all(FLERR,"Incorrect {} format in data file",keyword);
@ -692,7 +692,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
pitorsion_atom6[m][num_pitorsion[m]] = atom6;
num_pitorsion[m]++;
}
if ((m = atom->map(atom2)) >= 0) {
if (num_pitorsion[m] == PITORSIONMAX)
error->one(FLERR,"Too many PiTorsions for one atom");
@ -705,7 +705,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
pitorsion_atom6[m][num_pitorsion[m]] = atom6;
num_pitorsion[m]++;
}
if ((m = atom->map(atom3)) >= 0) {
if (num_pitorsion[m] == PITORSIONMAX)
error->one(FLERR,"Too many PiTorsions for one atom");
@ -731,7 +731,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
pitorsion_atom6[m][num_pitorsion[m]] = atom6;
num_pitorsion[m]++;
}
if ((m = atom->map(atom5)) >= 0) {
if (num_pitorsion[m] == PITORSIONMAX)
error->one(FLERR,"Too many PiTorsions for one atom");
@ -744,7 +744,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
pitorsion_atom6[m][num_pitorsion[m]] = atom6;
num_pitorsion[m]++;
}
if ((m = atom->map(atom6)) >= 0) {
if (num_pitorsion[m] == PITORSIONMAX)
error->one(FLERR,"Too many PiTorsions for one atom");
@ -757,7 +757,7 @@ void FixAmoebaPiTorsion::read_data_section(char *keyword, int n, char *buf,
pitorsion_atom6[m][num_pitorsion[m]] = atom6;
num_pitorsion[m]++;
}
buf = next + 1;
}
}
@ -780,7 +780,7 @@ bigint FixAmoebaPiTorsion::read_data_skip_lines(char *keyword)
void FixAmoebaPiTorsion::write_data_header(FILE *fp, int mth)
{
if (mth == 0) fprintf(fp,BIGINT_FORMAT " pitorsions\n",npitorsions);
else if (mth == 1)
else if (mth == 1)
fprintf(fp,BIGINT_FORMAT " pitorsion types\n",npitorsion_types);
}
@ -898,7 +898,7 @@ void FixAmoebaPiTorsion::write_data_section(int mth, FILE *fp,
if (mth == 0) {
for (int i = 0; i < n; i++)
fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
" " TAGINT_FORMAT "\n",
index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
(tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,

12
src/AMOEBA/improper_amoeba.cpp
View File

@ -127,9 +127,9 @@ void ImproperAmoeba::compute(int eflag, int vflag)
xcd = xic - xid;
ycd = yic - yid;
zcd = zic - zid;
// Allinger angle between A-C-D plane and D-B vector for D-B < AC
rad2 = xad*xad + yad*yad + zad*zad;
rcd2 = xcd*xcd + ycd*ycd + zcd*zcd;
dot = xad*xcd + yad*ycd + zad*zcd;
@ -152,12 +152,12 @@ void ImproperAmoeba::compute(int eflag, int vflag)
dt2 = dt * dt;
dt3 = dt2 * dt;
dt4 = dt2 * dt2;
e = eprefactor * k[type] * dt2 *
(1.0 + opbend_cubic*dt + opbend_quartic*dt2 +
e = eprefactor * k[type] * dt2 *
(1.0 + opbend_cubic*dt + opbend_quartic*dt2 +
opbend_pentic*dt3 + opbend_sextic*dt4);
deddt = fprefactor * k[type] * dt *
(2.0 + 3.0*opbend_cubic*dt + 4.0*opbend_quartic*dt2 +
deddt = fprefactor * k[type] * dt *
(2.0 + 3.0*opbend_cubic*dt + 4.0*opbend_quartic*dt2 +
5.0*opbend_pentic*dt3 + 6.0*opbend_sextic*dt4);
sign = (ee >= 0.0) ? 1.0 : -1.0;
dedcos = -deddt * sign / sqrt(cc*rdb2 - ee*ee);

383
src/AMOEBA/pair_amoeba.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -12,29 +12,30 @@
------------------------------------------------------------------------- */
#include "pair_amoeba.h"
#include <mpi.h>
#include "amoeba_convolution.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fft3d_wrap.h"
#include "fix.h"
#include "fix_store.h"
#include "force.h"
#include "gridcomm.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "update.h"
#include "utils.h"
#include <cstdlib>
#include <cstring>
#include <cctype>
#include "amoeba_convolution.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "fix_store.h"
#include "fft3d_wrap.h"
#include "gridcomm.h"
#include "my_page.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
using namespace LAMMPS_NS;
@ -56,7 +57,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
{
// error checks
if (strcmp(update->unit_style,"real") != 0)
if (strcmp(update->unit_style,"real") != 0)
error->all(FLERR,"Pair style amoeba/hippo require real units");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style amoeba/hippo require newton pair on");
@ -65,7 +66,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
me = comm->me;
nprocs = comm->nprocs;
// force field settings
one_coeff = 1;
@ -77,7 +78,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
optorder = 0;
maxualt = 0;
tcgnab = 0;
nmax = 0;
xaxis2local = yaxis2local = zaxis2local = NULL;
rpole = NULL;
@ -102,7 +103,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
m_kspace = p_kspace = pc_kspace = d_kspace = NULL;
i_kspace = ic_kspace = NULL;
numneigh_dipole = NULL;
firstneigh_dipole = NULL;
firstneigh_dipdip = NULL;
@ -135,7 +136,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
//electric = force->qqr2e;
// factors for FFT grid size
nfactors = 3;
factors = new int[nfactors];
factors[0] = 2;
@ -265,7 +266,7 @@ void PairAmoeba::compute(int eflag, int vflag)
if (atom->nlocal + atom->nghost > nmax) grow_local();
// set amtype/amgroup ptrs for rest of compute() to use it
amtype = atom->ivector[index_amtype];
amgroup = atom->ivector[index_amgroup];
@ -275,30 +276,30 @@ void PairAmoeba::compute(int eflag, int vflag)
// -------------------------------------------------------------------
// assignment of atoms to polarization groups
if (first_flag_compute) assign_groups();
// assigmment of multipole neighbors to each owned atom
// sets xaxis,yaxis,zaxis
// for HIPPO, also set pval for each atom, then ghost comm of pval
if (first_flag_compute) {
cfstyle = KMPOLE;
comm->forward_comm_pair(this);
comm->forward_comm(this);
kmpole();
if (hippo) {
pval = atom->dvector[index_pval];
double **pole = fixpole->astore;
int nlocal = atom->nlocal;
int itype,iclass;
for (int i = 0; i < nlocal; i++) {
itype = amtype[i];
iclass = amtype2class[itype];
pval[i] = pole[i][0] - pcore[iclass];
itype = amtype[i];
iclass = amtype2class[itype];
pval[i] = pole[i][0] - pcore[iclass];
}
cfstyle = PVAL;
comm->forward_comm_pair(this);
comm->forward_comm(this);
}
}
@ -328,7 +329,7 @@ void PairAmoeba::compute(int eflag, int vflag)
// reset KSpace recip matrix if box size/shape change dynamically
if (domain->box_change) lattice();
// compute reduced H coords for owned atoms
// needs to be computed before forward_comm with cfstyle = SETUP
@ -345,7 +346,7 @@ void PairAmoeba::compute(int eflag, int vflag)
rdn = kred[iclass];
xred[i][0] = rdn*(x[i][0]-x[j][0]) + x[j][0];
xred[i][1] = rdn*(x[i][1]-x[j][1]) + x[j][1];
xred[i][2] = rdn*(x[i][2]-x[j][2]) + x[j][2];
xred[i][2] = rdn*(x[i][2]-x[j][2]) + x[j][2];
}
}
@ -364,7 +365,7 @@ void PairAmoeba::compute(int eflag, int vflag)
if (amoeba) cfstyle = SETUP_AMOEBA;
else cfstyle = SETUP_HIPPO;
comm->forward_comm_pair(this);
comm->forward_comm(this);
if (amoeba) pbc_xred();
@ -400,7 +401,7 @@ void PairAmoeba::compute(int eflag, int vflag)
if (polar_rspace_flag || polar_kspace_flag) {
induce();
cfstyle = INDUCE;
comm->forward_comm_pair(this);
comm->forward_comm(this);
}
time6 = MPI_Wtime();
@ -424,14 +425,14 @@ void PairAmoeba::compute(int eflag, int vflag)
int nall = nlocal + atom->nghost;
for (int i = 0; i < 6; i++)
virial[i] = virhal[i] + virrepulse[i] + virdisp[i] +
virial[i] = virhal[i] + virrepulse[i] + virdisp[i] +
virpolar[i] + virmpole[i] + virqxfer[i];
// virial computation
// NOTE: how does this work for AMOEBA ?
// do all terms get summed this way, or only pairwise
// it is 2x what summed virial[6] above is
// if (vflag_fdotr) virial_fdotr_compute();
// timing information
@ -527,7 +528,7 @@ void PairAmoeba::allocate()
void PairAmoeba::settings(int narg, char **arg)
{
// turn on all FF components by default
hal_flag = repulse_flag = qxfer_flag = 1;
disp_rspace_flag = disp_kspace_flag = 1;
polar_rspace_flag = polar_kspace_flag = 1;
@ -565,15 +566,15 @@ void PairAmoeba::settings(int narg, char **arg)
if (strcmp(arg[iarg],"hal") == 0) hal_flag = newvalue;
else if (strcmp(arg[iarg],"repulse") == 0) repulse_flag = newvalue;
else if (strcmp(arg[iarg],"qxfer") == 0) qxfer_flag = newvalue;
else if (strcmp(arg[iarg],"disp") == 0)
else if (strcmp(arg[iarg],"disp") == 0)
disp_rspace_flag = disp_kspace_flag = newvalue;
else if (strcmp(arg[iarg],"disp/rspace") == 0) disp_rspace_flag = newvalue;
else if (strcmp(arg[iarg],"disp/kspace") == 0) disp_kspace_flag = newvalue;
else if (strcmp(arg[iarg],"polar") == 0)
else if (strcmp(arg[iarg],"polar") == 0)
polar_rspace_flag = polar_kspace_flag = newvalue;
else if (strcmp(arg[iarg],"polar/rspace") == 0) polar_rspace_flag = newvalue;
else if (strcmp(arg[iarg],"polar/kspace") == 0) polar_kspace_flag = newvalue;
else if (strcmp(arg[iarg],"mpole") == 0)
else if (strcmp(arg[iarg],"mpole") == 0)
mpole_rspace_flag = mpole_kspace_flag = newvalue;
else if (strcmp(arg[iarg],"mpole/rspace") == 0) mpole_rspace_flag = newvalue;
else if (strcmp(arg[iarg],"mpole/kspace") == 0) mpole_kspace_flag = newvalue;
@ -693,44 +694,44 @@ void PairAmoeba::init_style()
if (apewald == aeewald && polar_kspace_flag && !mpole_kspace_flag)
error->all(FLERR,
"Pair amoeba with apewald = aeewald requires mpole and polar together");
"Pair amoeba with apewald = aeewald requires mpole and polar together");
// check if all custom atom arrays were set via fix property/atom
int flag,cols;
index_amtype = atom->find_custom("amtype",flag,cols);
if (index_amtype < 0 || flag || cols)
if (index_amtype < 0 || flag || cols)
error->all(FLERR,"Pair amoeba amtype is not defined");
index_amgroup = atom->find_custom("amgroup",flag,cols);
if (index_amgroup < 0 || flag || cols)
if (index_amgroup < 0 || flag || cols)
error->all(FLERR,"Pair amoeba amgroup is not defined");
index_ired = atom->find_custom("ired",flag,cols);
if (index_ired < 0 || flag || cols)
if (index_ired < 0 || flag || cols)
error->all(FLERR,"Pair amoeba ired is not defined");
index_xaxis = atom->find_custom("xaxis",flag,cols);
if (index_xaxis < 0 || flag || cols)
if (index_xaxis < 0 || flag || cols)
error->all(FLERR,"Pair amoeba xaxis is not defined");
index_yaxis = atom->find_custom("yaxis",flag,cols);
if (index_yaxis < 0 || flag || cols)
if (index_yaxis < 0 || flag || cols)
error->all(FLERR,"Pair amoeba yaxis is not defined");
index_zaxis = atom->find_custom("zaxis",flag,cols);
if (index_zaxis < 0 || flag || cols)
if (index_zaxis < 0 || flag || cols)
error->all(FLERR,"Pair amoeba zaxis is not defined");
index_polaxe = atom->find_custom("polaxe",flag,cols);
if (index_polaxe < 0 || flag || cols)
if (index_polaxe < 0 || flag || cols)
error->all(FLERR,"Pair amoeba polaxe is not defined");
index_pval = atom->find_custom("pval",flag,cols);
if (index_pval < 0 || !flag || cols)
if (index_pval < 0 || !flag || cols)
error->all(FLERR,"Pair amoeba pval is not defined");
// -------------------------------------------------------------------
// one-time initializations
// can't do earlier b/c need all atoms to exist
// -------------------------------------------------------------------
// creation of per-atom storage
// create a new fix STORE style for each atom's pole vector
// id = "AMOEBA_pole", fix group = all
@ -815,24 +816,24 @@ void PairAmoeba::init_style()
// initialize KSpace Ewald settings and FFTs and parallel grid objects
// Coulombic grid is used with two orders: bseorder and bsporder
// so need two GridComm instantiations for ghost comm
if (first_flag) {
kewald();
if (use_ewald) {
m_kspace =
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bseorder,MPOLE_GRID);
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bseorder,MPOLE_GRID);
p_kspace =
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,POLAR_GRID);
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,POLAR_GRID);
pc_kspace =
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,POLAR_GRIDC);
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,POLAR_GRIDC);
i_kspace =
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,INDUCE_GRID);
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,INDUCE_GRID);
ic_kspace =
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,INDUCE_GRIDC);
new AmoebaConvolution(lmp,this,nefft1,nefft2,nefft3,bsporder,INDUCE_GRIDC);
}
if (use_dewald) {
d_kspace =
new AmoebaConvolution(lmp,this,ndfft1,ndfft2,ndfft3,bsdorder,DISP_GRID);
new AmoebaConvolution(lmp,this,ndfft1,ndfft2,ndfft3,bsdorder,DISP_GRID);
}
}
@ -861,7 +862,7 @@ void PairAmoeba::init_style()
csixpr = 0.0;
for (int i = 1; i <= n_amclass; i++) {
for (int j = i+1; j <= n_amclass; j++) {
csixpr += csix[i]*csix[j] * csix_num[i]*csix_num[j];
csixpr += csix[i]*csix[j] * csix_num[i]*csix_num[j];
}
}
csixpr *= 2.0;
@ -882,9 +883,9 @@ void PairAmoeba::init_style()
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) pval[i] = 0.0;
}
// output FF settings to screen and logfile
if (first_flag) {
if (comm->me == 0) {
if (screen) print_settings(screen);
@ -893,7 +894,7 @@ void PairAmoeba::init_style()
}
// all done with one-time initializations
first_flag = 0;
// -------------------------------------------------------------------
@ -927,7 +928,7 @@ void PairAmoeba::init_style()
ired = atom->ivector[index_ired];
int **nspecial = atom->nspecial;
int **special = atom->special;
tagint **special = atom->special;
int nlocal = atom->nlocal;
int itype,iclass;
@ -998,14 +999,14 @@ void PairAmoeba::print_settings(FILE *fp)
sqrt(off2),sqrt(cut2),
special_repel[1],special_repel[2],special_repel[3],special_repel[4]);
}
if (hippo) {
choose(QFER);
fprintf(fp," qxfer: cut %g taper %g mscale %g %g %g %g\n",
sqrt(off2),sqrt(cut2),
special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
}
if (hippo) {
if (use_dewald) {
choose(DISP_LONG);
@ -1020,47 +1021,47 @@ void PairAmoeba::print_settings(FILE *fp)
special_disp[1],special_disp[2],special_disp[3],special_disp[4]);
}
}
if (use_ewald) {
choose(MPOLE_LONG);
fprintf(fp," multipole: cut %g aewald %g bsorder %d "
"FFT %d %d %d mscale %g %g %g %g\n",
sqrt(off2),aewald,bseorder,nefft1,nefft2,nefft3,
special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
"FFT %d %d %d mscale %g %g %g %g\n",
sqrt(off2),aewald,bseorder,nefft1,nefft2,nefft3,
special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
} else {
choose(MPOLE);
fprintf(fp," multipole: cut %g aewald %g mscale %g %g %g %g\n",
sqrt(off2),aewald,
special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
sqrt(off2),aewald,
special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
}
if (use_ewald) {
choose(POLAR_LONG);
fprintf(fp," polar: cut %g aewald %g bsorder %d FFT %d %d %d\n",
sqrt(off2),aewald,bsporder,nefft1,nefft2,nefft3);
sqrt(off2),aewald,bsporder,nefft1,nefft2,nefft3);
fprintf(fp," pscale %g %g %g %g piscale %g %g %g %g "
"wscale %g %g %g %g d/u scale %g %g\n",
special_polar_pscale[1],special_polar_pscale[2],
special_polar_pscale[3],special_polar_pscale[4],
special_polar_piscale[1],special_polar_piscale[2],
special_polar_piscale[3],special_polar_piscale[4],
special_polar_wscale[1],special_polar_wscale[2],
special_polar_wscale[3],special_polar_wscale[4],
polar_dscale,polar_uscale);
"wscale %g %g %g %g d/u scale %g %g\n",
special_polar_pscale[1],special_polar_pscale[2],
special_polar_pscale[3],special_polar_pscale[4],
special_polar_piscale[1],special_polar_piscale[2],
special_polar_piscale[3],special_polar_piscale[4],
special_polar_wscale[1],special_polar_wscale[2],
special_polar_wscale[3],special_polar_wscale[4],
polar_dscale,polar_uscale);
} else {
choose(POLAR);
fprintf(fp," polar: cut %g aewald %g\n",sqrt(off2),aewald);
fprintf(fp," pscale %g %g %g %g piscale %g %g %g %g "
"wscale %g %g %g %g d/u scale %g %g\n",
special_polar_pscale[1],special_polar_pscale[2],
special_polar_pscale[3],special_polar_pscale[4],
special_polar_piscale[1],special_polar_piscale[2],
special_polar_piscale[3],special_polar_piscale[4],
special_polar_wscale[1],special_polar_wscale[2],
special_polar_wscale[3],special_polar_wscale[4],
polar_dscale,polar_uscale);
"wscale %g %g %g %g d/u scale %g %g\n",
special_polar_pscale[1],special_polar_pscale[2],
special_polar_pscale[3],special_polar_pscale[4],
special_polar_piscale[1],special_polar_piscale[2],
special_polar_piscale[3],special_polar_piscale[4],
special_polar_wscale[1],special_polar_wscale[2],
special_polar_wscale[3],special_polar_wscale[4],
polar_dscale,polar_uscale);
}
choose(USOLV);
fprintf(fp," precondition: cut %g\n",sqrt(off2));
}
@ -1191,7 +1192,7 @@ int PairAmoeba::pack_forward_comm(int n, int *list, double *buf,
} else if (cfstyle == KMPOLE) {
int **nspecial = atom->nspecial;
int **special = atom->special;
tagint **special = atom->special;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(nspecial[j][0]).d;
@ -1296,7 +1297,7 @@ void PairAmoeba::unpack_forward_comm(int n, int first, double *buf)
} else if (cfstyle == KMPOLE) {
int **nspecial = atom->nspecial;
int **special = atom->special;
tagint **special = atom->special;
for (i = first; i < last; i++) {
nspecial[i][0] = (int) ubuf(buf[m++]).i;
for (k = 0; k < nspecial[i][0]; k++)
@ -1598,7 +1599,7 @@ void PairAmoeba::assign_groups()
int maxstack = 0;
int *stack = NULL;
int **special = atom->special;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
@ -1606,91 +1607,91 @@ void PairAmoeba::assign_groups()
// initially, groupID = atomID
// communicate new groupIDs to ghost atoms
for (i = 0; i < nlocal; i++) amgroup[i] = tag[i];
cfstyle = AMGROUP;
comm->forward_comm_pair(this);
comm->forward_comm(this);
// loop until no ghost atom groupIDs are reset
while (1) {
// loop over all atoms and their group neighborhoods
ghostmark = 0;
for (i = 0; i < nlocal; i++) {
// push atom I on stack
nstack = 0;
if (nstack == maxstack) {
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
}
stack[nstack++] = i;
// loop over I's group neighborhood until stack is empty
while (nstack > 0) {
// pop atom M off stack
m = stack[nstack-1];
nstack--;
mtype = amtype[m];
// pop atom M off stack
// loop over bond partners of atom M
nbond = nspecial[m][0];
for (ibond = 0; ibond < nbond; ibond++) {
jglobal = special[m][ibond];
j = atom->map(jglobal);
if (j < 0)
error->one(FLERR,"AMOEBA group assignment bond neighbor not found");
jtype = amtype[j];
m = stack[nstack-1];
nstack--;
mtype = amtype[m];
// if amtype of bondpartner J is not in polgroup of atom M, continue
ngroup = npolgroup[mtype];
for (igroup = 0; igroup < ngroup; igroup++)
if (jtype == polgroup[mtype][igroup]) break;
if (igroup == ngroup) continue;
// if groupID of atoms J and M are the same, continue
// else set atom with larger groupID to smaller groupID
// if changed atom is ghost, set ghostmark, else push atom on stack
if (amgroup[m] == amgroup[j]) continue;
// loop over bond partners of atom M
if (amgroup[m] > amgroup[j]) {
amgroup[m] = amgroup[j];
if (nstack == maxstack) {
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
}
stack[nstack++] = m;
} else {
amgroup[j] = amgroup[m];
if (j >= nlocal) ghostmark = 1;
else {
if (nstack == maxstack) {
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
}
stack[nstack++] = j;
}
}
}
nbond = nspecial[m][0];
for (ibond = 0; ibond < nbond; ibond++) {
jglobal = special[m][ibond];
j = atom->map(jglobal);
if (j < 0)
error->one(FLERR,"AMOEBA group assignment bond neighbor not found");
jtype = amtype[j];
// if amtype of bondpartner J is not in polgroup of atom M, continue
ngroup = npolgroup[mtype];
for (igroup = 0; igroup < ngroup; igroup++)
if (jtype == polgroup[mtype][igroup]) break;
if (igroup == ngroup) continue;
// if groupID of atoms J and M are the same, continue
// else set atom with larger groupID to smaller groupID
// if changed atom is ghost, set ghostmark, else push atom on stack
if (amgroup[m] == amgroup[j]) continue;
if (amgroup[m] > amgroup[j]) {
amgroup[m] = amgroup[j];
if (nstack == maxstack) {
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
}
stack[nstack++] = m;
} else {
amgroup[j] = amgroup[m];
if (j >= nlocal) ghostmark = 1;
else {
if (nstack == maxstack) {
maxstack += DELTASTACK;
memory->grow(stack,maxstack,"amoeba:stack");
}
stack[nstack++] = j;
}
}
}
}
}
// communicate new groupIDs to ghost atoms
cfstyle = AMGROUP;
comm->forward_comm_pair(this);
comm->forward_comm(this);
// done if no proc reset groupID of a ghost atom
MPI_Allreduce(&ghostmark,&anyghostmark,1,MPI_INT,MPI_MAX,world);
if (!anyghostmark) break;
}
@ -1721,7 +1722,7 @@ void PairAmoeba::assign_groups()
void PairAmoeba::pbc_xred()
{
double prd,prd_half,delta;
double **x = atom->x;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
@ -1732,9 +1733,9 @@ void PairAmoeba::pbc_xred()
for (int i = nlocal; i < nall; i++) {
delta = xred[i][0] - x[i][0];
while (fabs(delta) > prd_half) {
if (delta < 0.0) xred[i][0] += prd;
else xred[i][0] -= prd;
delta = xred[i][0] - x[i][0];
if (delta < 0.0) xred[i][0] += prd;
else xred[i][0] -= prd;
delta = xred[i][0] - x[i][0];
}
}
}
@ -1745,22 +1746,22 @@ void PairAmoeba::pbc_xred()
for (int i = nlocal; i < nall; i++) {
delta = xred[i][1] - x[i][1];
while (fabs(delta) > prd_half) {
if (delta < 0.0) xred[i][1] += prd;
else xred[i][1] -= prd;
delta = xred[i][1] - x[i][1];
if (delta < 0.0) xred[i][1] += prd;
else xred[i][1] -= prd;
delta = xred[i][1] - x[i][1];
}
}
}
if (domain->zperiodic) {
prd = domain->zprd;
prd_half = domain->zprd_half;
for (int i = nlocal; i < nall; i++) {
delta = xred[i][2] - x[i][2];
while (fabs(delta) > prd_half) {
if (delta < 0.0) xred[i][2] += prd;
else xred[i][2] -= prd;
delta = xred[i][2] - x[i][2];
if (delta < 0.0) xred[i][2] += prd;
else xred[i][2] -= prd;
delta = xred[i][2] - x[i][2];
}
}
}
@ -1846,7 +1847,7 @@ void PairAmoeba::allocate_vdwl()
memory->create(epsilon,n_amclass+1,n_amclass+1,"amoeba:epsilon");
memory->create(epsilon4,n_amclass+1,n_amclass+1,"amoeba:epsilon4");
}
void PairAmoeba::deallocate_vdwl()
{
memory->destroy(radmin);
@ -1922,7 +1923,7 @@ void PairAmoeba::choose(int which)
double off,cut;
// short-range only terms
if (which == HAL) {
off = vdwcut;
cut = vdwtaper;
@ -1949,7 +1950,7 @@ void PairAmoeba::choose(int which)
cut = 0.99*off; // not used
// short-range + long-range terms
} else if (which == DISP_LONG) {
off = dispcut;
cut = 0.99*cut; // not used
@ -1963,7 +1964,7 @@ void PairAmoeba::choose(int which)
cut = 0.99*off; // not used
aewald = apewald;
}
off2 = off*off;
cut2 = cut*cut;
@ -1993,16 +1994,16 @@ void PairAmoeba::mix()
double TWOSIX = pow(2.0,1.0/6.0);
for (i = 1; i <= n_amclass; i++) {
// printf("MIX i %d nclass %d eps %g sigma %g\n",
// i,n_amclass,vdwl_eps[i],vdwl_sigma[i]);
for (j = i; j <= n_amclass; j++) {
ei = vdwl_eps[i];
ej = vdwl_eps[j];
ri = vdwl_sigma[i];
rj = vdwl_sigma[j];
if (radius_type == SIGMA) {
ri *= TWOSIX;
rj *= TWOSIX;
@ -2011,7 +2012,7 @@ void PairAmoeba::mix()
ri *= 0.5;
rj *= 0.5;
}
sri = sqrt(ri);
ei = fabs(ei);
sei = sqrt(ei);
@ -2061,7 +2062,7 @@ void PairAmoeba::mix()
j = vdwl_class_pair[m][1];
rij = vdwl_sigma_pair[m];
eij = vdwl_eps_pair[m];
if (radius_type == SIGMA) rij *= TWOSIX;
radmin[j][i] = radmin[i][j] = rij;
@ -2217,12 +2218,12 @@ void PairAmoeba::grow_local()
------------------------------------------------------------------------- */
void PairAmoeba::dump6(FILE *fp, const char *columns, double scale,
double **a, double **b)
double **a, double **b)
{
int i,j,m;
MPI_Status status;
MPI_Request request;
// setup
int size_one = 7;
@ -2253,7 +2254,7 @@ void PairAmoeba::dump6(FILE *fp, const char *columns, double scale,
}
// write file
if (me == 0) {
fprintf(fp,"ITEM: TIMESTEP\n");
fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep);
@ -2268,25 +2269,25 @@ void PairAmoeba::dump6(FILE *fp, const char *columns, double scale,
int nlines;
double tmp;
if (me == 0) {
for (int iproc = 0; iproc < nprocs; iproc++) {
if (iproc) {
MPI_Irecv(buf,maxlocal*size_one,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= size_one;
MPI_Irecv(buf,maxlocal*size_one,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= size_one;
} else nlines = nlocal;
m = 0;
for (i = 0; i < nlines; i++) {
for (j = 0; j < size_one; j++) {
if (j == 0) fprintf(fp,"%d",static_cast<tagint> (buf[m]));
else fprintf(fp," %g",buf[m]);
m++;
}
fprintf(fp,"\n");
for (j = 0; j < size_one; j++) {
if (j == 0) fprintf(fp,"%d",static_cast<tagint> (buf[m]));
else fprintf(fp," %g",buf[m]);
m++;
}
fprintf(fp,"\n");
}
}

50
src/AMOEBA/pair_amoeba.h
View File

@ -1,6 +1,6 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -86,7 +86,7 @@ class PairAmoeba : public Pair {
double time_init,time_hal,time_repulse,time_disp;
double time_mpole,time_induce,time_polar,time_qxfer;
// energy/virial components
double ehal,erepulse,edisp,epolar,empole,eqxfer;
@ -116,7 +116,7 @@ class PairAmoeba : public Pair {
double special_polar_wscale[5];
double polar_dscale,polar_uscale;
// scalar values defined in keyfile
double dhal,ghal;
@ -132,7 +132,7 @@ class PairAmoeba : public Pair {
double usolvcut;
int use_ewald,use_dewald;
int use_pred;
int politer,polpred;
int pcgprec,pcgguess;
@ -144,7 +144,7 @@ class PairAmoeba : public Pair {
int aeewald_key,apewald_key,adewald_key;
int pmegrid_key,dpmegrid_key;
// types and classes
int n_amtype; // # of defined AMOEBA types, 1-N
@ -155,7 +155,7 @@ class PairAmoeba : public Pair {
int *amtype_defined; // 1 if type was defined in FF file
int *amclass_defined; // 1 if class was defined in FF file
int *amtype2class; // amt2c[i] = class which type I belongs to
// static per-atom properties, must persist as atoms migrate
int index_amtype,index_amgroup,index_ired;
@ -168,7 +168,7 @@ class PairAmoeba : public Pair {
tagint *xaxis,*yaxis,*zaxis; // IDs of nearby atoms for multipole def
int *polaxe;
double *pval;
char *id_pole,*id_udalt,*id_upalt;
class FixStore *fixpole; // stores pole = multipole components
class FixStore *fixudalt; // stores udalt = induced dipole history
@ -194,7 +194,7 @@ class PairAmoeba : public Pair {
// multipole frame info for each amtype, read from PRM file
int *nmultiframe; // # of frames for each type
int **mpaxis; // polaxe values
int **mpaxis; // polaxe values
int **xpole,**ypole,**zpole; // other types in xyz dirs for multipole frame
double ***fpole; // 13 values from file
// 0 = monopole, same as q
@ -216,7 +216,7 @@ class PairAmoeba : public Pair {
double *vdwl_sigma_pair; // Vdwl sigma for pair of classes
int nvdwl_pair; // # of pairwise Vdwl entries in file
int max_vdwl_pair; // size of allocated data for pairwise Vdwl
// vectors and arrays of small size
double *copt,*copm; // 0:optorder in length
@ -245,7 +245,7 @@ class PairAmoeba : public Pair {
// peratom values computed each step
// none of them persist with atoms
// some of them need communication to ghosts
double **rpole; // multipole, comm to ghosts
int *xaxis2local,*yaxis2local,*zaxis2local; // xyz axis IDs -> local indices
@ -271,7 +271,7 @@ class PairAmoeba : public Pair {
double ***fopt,***foptp; // computed in induce, used by polar, if OPT
// Nlocal x optorder x 10
// derived local neighbor lists
int *numneigh_dipole; // number of dipole neighs for each atom
@ -292,7 +292,7 @@ class PairAmoeba : public Pair {
// in indices = owned portion of grid in spatial decomp
// out indices = in + ghost grid cells
// fft indices = owned portion of grid in FFT decomp
int nefft1,nefft2,nefft3; // for electrostatic PME operations
int ndfft1,ndfft2,ndfft3; // for dispersion PME operations
@ -300,12 +300,12 @@ class PairAmoeba : public Pair {
int bsporder; // for polarization
int bsdorder; // for dispersion
int bsordermax; // max of 3 bsorder values
double aewald; // current Ewald alpha
double aeewald; // for electrostatics
double apewald; // for polarization
double adewald; // for dispersion
double *bsmod1,*bsmod2,*bsmod3; // B-spline module along abc axes
// set to max of any nfft1,nfft2,nfft3
@ -318,12 +318,12 @@ class PairAmoeba : public Pair {
// indices ARE flipped vs Fortran
// Kspace data for induce and polar
double *qfac; // convoulution pre-factors
double **cmp,**fmp,**cphi,**fphi; // Cartesian and fractional multipoles
double **cmp,**fmp,**cphi,**fphi; // Cartesian and fractional multipoles
// params for current KSpace solve and FFT being worked on
int nfft1,nfft2,nfft3; // size of FFT
int bsorder; // stencil size
double recip[3][3]; // indices NOT flipped vs Fortran
@ -332,12 +332,12 @@ class PairAmoeba : public Pair {
class AmoebaConvolution *m_kspace,*p_kspace,*pc_kspace,*d_kspace;
class AmoebaConvolution *i_kspace,*ic_kspace;
// FFT grid size factors
int nfactors; // # of factors
int *factors; // list of possible factors (2,3,5)
// components of force field
void hal();
@ -396,9 +396,9 @@ class PairAmoeba : public Pair {
void kewald();
void kewald_parallel(int, int, int, int,
int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &);
int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &);
double ewaldcof(double);
int factorable(int);
@ -411,7 +411,7 @@ class PairAmoeba : public Pair {
void allocate();
void print_settings(FILE *);
void initialize_vdwl();
void allocate_vdwl();
void deallocate_vdwl();
@ -447,7 +447,7 @@ class PairAmoeba : public Pair {
void set_defaults();
void read_prmfile(char *);
void read_keyfile(char *);
int read_section_name(FILE *fp, char *);
int read_section_line(FILE *fp, char *, int &, char *);
int tokenize(char *, char **&, char *&);
@ -477,7 +477,7 @@ class PairAmoeba : public Pair {
void file_charge_transfer(int, char **);
// inline function for neighbor list unmasking
inline int sbmask15(int j) const {
return j >> SBBITS15 & 7;
}

2
src/AMOEBA/pair_hippo.cpp
View File

@ -1,6 +1,6 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract

2
src/AMOEBA/pair_hippo.h
View File

@ -1,6 +1,6 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract

1
src/Makefile
View File

@ -146,6 +146,7 @@ PACKAGE = \
PACKBASIC = kspace manybody molecule rigid
PACKMOST = \
amoeba \
asphere \
bocs \
body \

2
src/angle.cpp
View File

@ -357,7 +357,7 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond, double ea
------------------------------------------------------------------------- */
void Angle::ev_tally4(int i, int j, int k, int m, int nlocal, int newton_bond,
double eangle,
double eangle,
double *f1, double *f2, double *f3, double *f4)
{
double eanglefourth,v[6];

4
src/atom.cpp
View File

@ -215,7 +215,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
maxspecial15 = 1;
nspecial15 = nullptr;
special15 = nullptr;
// DIELECTRIC package
area = ed = em = epsilon = curvature = q_unscaled = nullptr;
@ -679,7 +679,7 @@ void Atom::set_atomflag_defaults()
contact_radius_flag = smd_data_9_flag = smd_stress_flag = 0;
eff_plastic_strain_flag = eff_plastic_strain_rate_flag = 0;
nspecial15_flag = 0;
pdscale = 1.0;
}

2
src/atom.h
View File

@ -208,7 +208,7 @@ class Atom : protected Pointers {
// AMOEBA package
int nspecial15_flag;
// Peridynamics scale factor, used by dump cfg
double pdscale;

8
src/fix_property_atom.cpp
View File

@ -304,7 +304,7 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf, tagint
atom->iarray[index[j]][m][k] = values.next_int();
} else if (styles[j] == DARRAY) {
ncol = cols[j];
for (k = 0; k < ncol; k++)
for (k = 0; k < ncol; k++)
atom->darray[index[j]][m][k] = values.next_double();
}
}
@ -570,7 +570,7 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
}
}
@ -776,11 +776,11 @@ int FixPropertyAtom::pack_restart(int i, double *buf)
buf[m++] = atom->dvector[index[nv]][i];
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
for (k = 0; k < ncol; k++)
buf[m++] = atom->darray[index[nv]][i][k];
}
}

4
src/fix_store.cpp
View File

@ -105,7 +105,7 @@ vstore(nullptr), astore(nullptr), rbuf(nullptr)
if (flavor == GLOBAL) {
if (vecflag) {
for (int i = 0; i < n1; i++)
for (int i = 0; i < n1; i++)
vstore[i] = 0.0;
} else if (arrayflag) {
for (int i = 0; i < n1; i++)
@ -117,7 +117,7 @@ vstore(nullptr), astore(nullptr), rbuf(nullptr)
if (flavor == PERATOM) {
int nlocal = atom->nlocal;
if (vecflag) {
for (int i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
vstore[i] = 0.0;
} else if (arrayflag) {
for (int i = 0; i < nlocal; i++)

8
src/force.cpp
View File

@ -55,7 +55,7 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp)
special_angle = special_dihedral = 0;
special_onefive = 0;
special_extra = 0;
dielectric = 1.0;
qqr2e_lammps_real = 332.06371; // these constants are toggled
qqr2e_charmm_real = 332.0716; // by new CHARMM pair styles
@ -706,7 +706,7 @@ void Force::set_special(int narg, char **arg)
special_coul[1] = special_coul[2] = special_coul[3] = 0.0;
special_angle = special_dihedral = 0;
special_onefive = 0;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "amber") == 0) {
@ -777,8 +777,8 @@ void Force::set_special(int narg, char **arg)
else if (strcmp(arg[iarg+1],"yes") == 0) special_onefive = 1;
else error->all(FLERR,"Illegal special_bonds command");
if (special_onefive && atom->nspecial15_flag == 0)
error->all(FLERR,"Cannot set special_bonds one/five if "
"atom style does not support it");
error->all(FLERR,"Cannot set special_bonds one/five if "
"atom style does not support it");
iarg += 2;
} else error->all(FLERR,"Illegal special_bonds command");
}

2
src/force.h
View File

@ -116,7 +116,7 @@ class Force : protected Pointers {
// 1 if only weight 1,4 atoms if in a dihedral
int special_extra; // extra space for added bonds
int special_onefive; // 0 if 1-5 neighbors are not stored, 1 if yes
Force(class LAMMPS *);
~Force() override;
void init();

14
src/memory.h
View File

@ -649,8 +649,8 @@ class Memory : protected Pointers {
template <typename TYPE>
TYPE ****create4d_offset_last(TYPE ****&array, int n1lo, int n1hi,
int n2lo, int n2hi, int n3lo, int n3hi, int n4,
const char *name)
int n2lo, int n2hi, int n3lo, int n3hi, int n4,
const char *name)
{
if (n1lo > n1hi || n2lo > n2hi || n3lo > n3hi || n4 <= 0) return nullptr;
@ -668,9 +668,9 @@ class Memory : protected Pointers {
template <typename TYPE>
TYPE ****create4d_offset_last(TYPE *****& /*array*/, int /*n1lo*/, int /*n1hi*/,
int /*n2lo*/, int /*n2hi*/,
int /*n3lo*/, int /*n3hi*/, int /*n4*/,
const char *name)
int /*n2lo*/, int /*n2hi*/,
int /*n3lo*/, int /*n3hi*/, int /*n4*/,
const char *name)
{fail(name); return nullptr;}
/* ----------------------------------------------------------------------
@ -679,7 +679,7 @@ class Memory : protected Pointers {
template <typename TYPE>
void destroy4d_offset_last(TYPE ****&array,
int n1_offset, int n2_offset, int n3_offset)
int n1_offset, int n2_offset, int n3_offset)
{
if (array == nullptr) return;
sfree(&array[n1_offset][n2_offset][n3_offset][0]);
@ -688,7 +688,7 @@ class Memory : protected Pointers {
sfree(&array[n1_offset]);
array = nullptr;
}
/* ----------------------------------------------------------------------
create a 5d array
------------------------------------------------------------------------- */

2
src/modify.cpp
View File

@ -783,7 +783,7 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
const char *exceptions[] =
{"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx",
"deprecated", "STORE/KIM", "amoeba/pitorsion", "amoeba/bitorsion",
"deprecated", "STORE/KIM", "amoeba/pitorsion", "amoeba/bitorsion",
nullptr};
if (domain->box_exist == 0) {

4
src/pair.h
View File

@ -192,12 +192,12 @@ class Pair : protected Pointers {
virtual void unpack_forward_comm(int, int, double *) {}
virtual int pack_reverse_comm(int, int, double *) { return 0; }
virtual void unpack_reverse_comm(int, int *, double *) {}
virtual void pack_forward_grid(int, void *, int, int *) {}
virtual void unpack_forward_grid(int, void *, int, int *) {}
virtual void pack_reverse_grid(int, void *, int, int *) {}
virtual void unpack_reverse_grid(int, void *, int, int *) {}
virtual double memory_usage();
void set_copymode(int value) { copymode = value; }

44
src/special.cpp
View File

@ -73,7 +73,7 @@ void Special::build()
// set onefive_flag if special_bonds command set it
onefive_flag = force->special_onefive;
// initialize nspecial counters to 0
int **nspecial = atom->nspecial;
@ -89,7 +89,7 @@ void Special::build()
if (onefive_flag) {
for (int i = 0; i < nlocal; i++) nspecial15[i] = 0;
}
// setup atomIDs and procowner vectors in rendezvous decomposition
atom_owners();
@ -107,7 +107,7 @@ void Special::build()
// done if special_bond weights for 1-3, 1-4 are set to 1.0
// onefive_flag must also be off, else 1-4 is needed to create 1-5
if (!onefive_flag &&
force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0 &&
force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) {
@ -132,7 +132,7 @@ void Special::build()
// done if special_bond weights for 1-4 are set to 1.0
// onefive_flag must also be off, else 1-4 is needed to create 1-5
if (!onefive_flag &&
force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) {
dedup();
@ -671,7 +671,7 @@ void Special::dedup()
int unique;
// dedup onetwo
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
@ -703,7 +703,7 @@ void Special::dedup()
atom->map_one(tag[i],-1);
for (j = 0; j < unique; j++) atom->map_one(onethree[i][j],-1);
}
// dedup onefour
for (i = 0; i < nlocal; i++) {
@ -728,11 +728,11 @@ void Special::dedup()
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial15[i]; j++) {
m = onefive[i][j];
if (atom->map(m) < 0) {
onefive[i][unique++] = m;
atom->map_one(m,0);
}
m = onefive[i][j];
if (atom->map(m) < 0) {
onefive[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial15[i] = unique;
atom->map_one(tag[i],-1);
@ -784,7 +784,7 @@ void Special::combine()
int unique,unique15;
int maxspecial = 0;
int maxspecial15 = 0;
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
@ -816,11 +816,11 @@ void Special::combine()
if (onefive_flag) {
unique15 = 0;
for (j = 0; j < nspecial15[i]; j++) {
m = onefive[i][j];
if (atom->map(m) < 0) {
unique15++;
atom->map_one(m,0);
}
m = onefive[i][j];
if (atom->map(m) < 0) {
unique15++;
atom->map_one(m,0);
}
}
maxspecial15 = MAX(maxspecial15,unique15);
}
@ -927,11 +927,11 @@ void Special::combine()
if (onefive_flag) {
unique15 = 0;
for (j = 0; j < nspecial15[i]; j++) {
m = onefive[i][j];
if (atom->map(m) < 0) {
special15[i][unique15++] = m;
atom->map_one(m,0);
}
m = onefive[i][j];
if (atom->map(m) < 0) {
special15[i][unique15++] = m;
atom->map_one(m,0);
}
}
nspecial15[i] = unique15;
}