forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7513 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
9a9ed90766
commit
8c3a07c169
|
|
@ -4,6 +4,7 @@ if (test $1 = 1) then
|
|||
|
||||
cp angle_cosine_shift.cpp ..
|
||||
cp angle_cosine_shift_exp.cpp ..
|
||||
cp angle_dipole.cpp ..
|
||||
cp bond_harmonic_shift.cpp ..
|
||||
cp bond_harmonic_shift_cut.cpp ..
|
||||
cp compute_ackland_atom.cpp ..
|
||||
|
|
@ -21,6 +22,7 @@ if (test $1 = 1) then
|
|||
|
||||
cp angle_cosine_shift.h ..
|
||||
cp angle_cosine_shift_exp.h ..
|
||||
cp angle_dipole.h ..
|
||||
cp bond_harmonic_shift.h ..
|
||||
cp bond_harmonic_shift_cut.h ..
|
||||
cp compute_ackland_atom.h ..
|
||||
|
|
@ -40,6 +42,7 @@ elif (test $1 = 0) then
|
|||
|
||||
rm -f ../angle_cosine_shift.cpp
|
||||
rm -f ../angle_cosine_shift_exp.cpp
|
||||
rm -f ../angle_dipole.cpp
|
||||
rm -f ../bond_harmonic_shift.cpp
|
||||
rm -f ../bond_harmonic_shift_cut.cpp
|
||||
rm -f ../compute_ackland_atom.cpp
|
||||
|
|
@ -57,6 +60,7 @@ elif (test $1 = 0) then
|
|||
|
||||
rm -f ../angle_cosine_shift.h
|
||||
rm -f ../angle_cosine_shift_exp.h
|
||||
rm -f ../angle_dipole.h
|
||||
rm -f ../bond_harmonic_shift.h
|
||||
rm -f ../bond_harmonic_shift_cut.h
|
||||
rm -f ../compute_ackland_atom.h
|
||||
|
|
|
|||
|
|
@ -20,6 +20,7 @@ about the feature or its coding.
|
|||
|
||||
angle_style cosine/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
angle_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
angle_style dipole, Mario Orsi, orsimario at gmail.com, 10 Jan 12
|
||||
bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
|
||||
|
|
|
|||
|
|
@ -0,0 +1,207 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Mario Orsi (U Southampton), orsimario@gmail.com
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "angle_dipole.h"
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleDipole::AngleDipole(LAMMPS *lmp) : Angle(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleDipole::~AngleDipole()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(k);
|
||||
memory->destroy(gamma0);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleDipole::compute(int eflag, int vflag)
|
||||
{
|
||||
int iRef,iDip,iDummy,n,type;
|
||||
double delx,dely,delz;
|
||||
double eangle,tangle,f1[3],f3[3];
|
||||
double r,dr,cosGamma,deltaGamma,kdg,rmu;
|
||||
|
||||
eangle = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double **x = atom->x; // position vector
|
||||
double **mu = atom->mu; // point-dipole components and moment magnitude
|
||||
double **torque = atom->torque;
|
||||
int **anglelist = neighbor->anglelist;
|
||||
int nanglelist = neighbor->nanglelist;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
if (!newton_bond)
|
||||
error->all(FLERR,"'newton' flag for bonded interactions must be 'on'");
|
||||
|
||||
for (n = 0; n < nanglelist; n++) {
|
||||
iDip = anglelist[n][0]; // dipole whose orientation is to be restrained
|
||||
iRef = anglelist[n][1]; // reference atom toward which dipole will point
|
||||
iDummy = anglelist[n][2]; // dummy atom - irrelevant to the interaction
|
||||
type = anglelist[n][3];
|
||||
|
||||
delx = x[iRef][0] - x[iDip][0];
|
||||
dely = x[iRef][1] - x[iDip][1];
|
||||
delz = x[iRef][2] - x[iDip][2];
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
|
||||
r = sqrt(delx*delx + dely*dely + delz*delz);
|
||||
|
||||
rmu = r * mu[iDip][3];
|
||||
cosGamma = (mu[iDip][0]*delx+mu[iDip][1]*dely+mu[iDip][2]*delz) / rmu;
|
||||
deltaGamma = cosGamma - cos(gamma0[type]);
|
||||
kdg = k[type] * deltaGamma;
|
||||
|
||||
if (eflag) eangle = kdg * deltaGamma; // energy
|
||||
|
||||
tangle = 2.0 * kdg / rmu;
|
||||
|
||||
torque[iDip][0] += tangle * (dely*mu[iDip][2] - delz*mu[iDip][1]);
|
||||
torque[iDip][1] += tangle * (delz*mu[iDip][0] - delx*mu[iDip][2]);
|
||||
torque[iDip][2] += tangle * (delx*mu[iDip][1] - dely*mu[iDip][0]);
|
||||
|
||||
f1[0] = f1[1] = f1[2] = f3[0] = f3[1] = f3[2] = 0.0;
|
||||
|
||||
if (evflag) // tally energy (virial=0 because force=0)
|
||||
ev_tally(iRef,iDip,iDummy,nlocal,newton_bond,eangle,f1,f3,
|
||||
0.0,0.0,0.0,0.0,0.0,0.0);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleDipole::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(gamma0,n+1,"angle:gamma0");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more types
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleDipole::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
force->bounds(arg[0],atom->nangletypes,ilo,ihi);
|
||||
|
||||
double k_one = force->numeric(arg[1]);
|
||||
double gamma0_one = force->numeric(arg[2]);
|
||||
|
||||
// convert gamma0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = k_one;
|
||||
gamma0[i] = gamma0_one/180.0 * MY_PI;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
used by SHAKE
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double AngleDipole::equilibrium_angle(int i)
|
||||
{
|
||||
return gamma0[i];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleDipole::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&gamma0[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleDipole::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&gamma0[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&gamma0[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
used by ComputeAngleLocal
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double AngleDipole::single(int type, int iRef, int iDip, int iDummy)
|
||||
{
|
||||
double **x = atom->x; // position vector
|
||||
double **mu = atom->mu; // point-dipole components and moment magnitude
|
||||
|
||||
double delx = x[iRef][0] - x[iDip][0];
|
||||
double dely = x[iRef][1] - x[iDip][1];
|
||||
double delz = x[iRef][2] - x[iDip][2];
|
||||
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
|
||||
double r = sqrt(delx*delx + dely*dely + delz*delz);
|
||||
double rmu = r * mu[iDip][3];
|
||||
double cosGamma = (mu[iDip][0]*delx+mu[iDip][1]*dely+mu[iDip][2]*delz) / rmu;
|
||||
double deltaGamma = cosGamma - cos(gamma0[type]);
|
||||
double kdg = k[type] * deltaGamma;
|
||||
|
||||
return kdg * deltaGamma; // energy
|
||||
}
|
||||
|
|
@ -0,0 +1,56 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef ANGLE_CLASS
|
||||
|
||||
AngleStyle(dipole,AngleDipole)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_ANGLE_DIPOLE_H
|
||||
#define LMP_ANGLE_DIPOLE_H
|
||||
|
||||
#include "stdio.h"
|
||||
#include "angle.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class AngleDipole : public Angle {
|
||||
public:
|
||||
AngleDipole(class LAMMPS *);
|
||||
virtual ~AngleDipole();
|
||||
virtual void compute(int, int);
|
||||
void coeff(int, char **);
|
||||
double equilibrium_angle(int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
double single(int, int, int, int);
|
||||
|
||||
protected:
|
||||
double *k,*gamma0;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args for angle coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,98 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
angle_style dipole command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
angle_style dipole :pre
|
||||
|
||||
[Examples:]
|
||||
|
||||
angle_style dipole
|
||||
angle_coeff 6 2.1 180.0 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
The {dipole} angle style is used to control the orientation of a dipolar
|
||||
atom within a molecule "(Orsi)"_#Orsi. Specifically, the {dipole} angle
|
||||
style restrains the orientation of a point dipole mu_j (embedded in atom
|
||||
'j') with respect to a reference (bond) vector r_ij = r_i - r_j, where 'i'
|
||||
is another atom of the same molecule (typically, 'i' and 'j' are also
|
||||
covalently bonded).
|
||||
|
||||
It is convenient to define an angle gamma between the 'free' vector mu_j
|
||||
and the reference (bond) vector r_ij:
|
||||
|
||||
:c,image(Eqs/angle_dipole_gamma.jpg)
|
||||
|
||||
The {dipole} angle style uses the potential:
|
||||
|
||||
:c,image(Eqs/angle_dipole_potential.jpg)
|
||||
|
||||
where K is a rigidity constant and gamma0 is an equilibrium (reference)
|
||||
angle.
|
||||
|
||||
The torque on the dipole can be obtained by differentiating the
|
||||
potential using the 'chain rule' as in appendix C.3 of
|
||||
"(Allen)"_#Allen:
|
||||
|
||||
:c,image(Eqs/angle_dipole_torque.jpg)
|
||||
|
||||
Example: if gamma0 is set to 0 degrees, the torque generated by
|
||||
the potential will tend to align the dipole along the reference
|
||||
direction defined by the (bond) vector r_ij (in other words, mu_j is
|
||||
restrained to point towards atom 'i').
|
||||
|
||||
Note that the angle dipole potential does not give rise to any force,
|
||||
because it does not depend on the distance between i and j (it only
|
||||
depends on the angle between mu_j and r_ij).
|
||||
|
||||
The following coefficients must be defined for each angle type via the
|
||||
"angle_coeff"_angle_coeff.html command as in the example above, or in
|
||||
the data file or restart files read by the "read_data"_read_data.html
|
||||
or "read_restart"_read_restart.html commands:
|
||||
|
||||
K (energy)
|
||||
gamma0 (degrees) :ul
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This angle style can only be used if LAMMPS was built with the
|
||||
"user-misc" package. See the "Making LAMMPS"_Section_start.html#2_3
|
||||
section for more info on packages.
|
||||
|
||||
IMPORTANT: In the "Angles" section of the data file, the atom id 'j'
|
||||
corresponding to the dipole to restrain must come before the atom id
|
||||
of the reference atom 'i'. A third atom id 'k' must also be provided,
|
||||
although 'k' is just a 'dummy' atom which can be any atom; it may be
|
||||
useful to choose a convention (e.g., 'k'='i') and adhere to it. For
|
||||
example, if id=1 for the dipolar atom to restrain, and id=2 for the
|
||||
reference atom, the corresponding line in the "Angles" section of the
|
||||
data file would read: X X 1 2 2
|
||||
|
||||
The "newton" command for intramolecular interactions must be "on"
|
||||
(which is the default).
|
||||
|
||||
This angle style should not be used with SHAKE.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"angle_coeff"_angle_coeff.html, "angle_hybrid"_angle_hybrid.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Orsi)
|
||||
[(Orsi)] Orsi & Essex, The ELBA force field for coarse-grain modeling of
|
||||
lipid membranes, PloS ONE 6(12): e28637, 2011.
|
||||
|
||||
:link(Allen)
|
||||
[(Allen)] Allen & Tildesley, Computer Simulation of Liquids,
|
||||
Clarendon Press, Oxford, 1987.
|
||||
Loading…
Reference in New Issue