forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5122 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
parent
4f64b5c792
commit
d6cfefc04c
|
@ -5,6 +5,7 @@ if (test $1 = 1) then
|
|||
cp angle_charmm.cpp ..
|
||||
cp angle_cosine.cpp ..
|
||||
cp angle_cosine_delta.cpp ..
|
||||
cp angle_cosine_periodic.cpp ..
|
||||
cp angle_cosine_squared.cpp ..
|
||||
cp angle_harmonic.cpp ..
|
||||
cp angle_hybrid.cpp ..
|
||||
|
@ -32,12 +33,16 @@ if (test $1 = 1) then
|
|||
cp improper_cvff.cpp ..
|
||||
cp improper_harmonic.cpp ..
|
||||
cp improper_hybrid.cpp ..
|
||||
cp improper_umbrella.cpp ..
|
||||
cp pair_hbond_dreiding_lj.cpp ..
|
||||
cp pair_hbond_dreiding_morse.cpp ..
|
||||
cp pair_lj_charmm_coul_charmm.cpp ..
|
||||
cp pair_lj_charmm_coul_charmm_implicit.cpp ..
|
||||
|
||||
cp angle_charmm.h ..
|
||||
cp angle_cosine.h ..
|
||||
cp angle_cosine_delta.h ..
|
||||
cp angle_cosine_periodic.h ..
|
||||
cp angle_cosine_squared.h ..
|
||||
cp angle_harmonic.h ..
|
||||
cp angle_hybrid.h ..
|
||||
|
@ -65,6 +70,9 @@ if (test $1 = 1) then
|
|||
cp improper_cvff.h ..
|
||||
cp improper_harmonic.h ..
|
||||
cp improper_hybrid.h ..
|
||||
cp improper_umbrella.h ..
|
||||
cp pair_hbond_dreiding_lj.h ..
|
||||
cp pair_hbond_dreiding_morse.h ..
|
||||
cp pair_lj_charmm_coul_charmm.h ..
|
||||
cp pair_lj_charmm_coul_charmm_implicit.h ..
|
||||
|
||||
|
@ -73,6 +81,7 @@ elif (test $1 = 0) then
|
|||
rm ../angle_charmm.cpp
|
||||
rm ../angle_cosine.cpp
|
||||
rm ../angle_cosine_delta.cpp
|
||||
rm ../angle_cosine_periodic.cpp
|
||||
rm ../angle_cosine_squared.cpp
|
||||
rm ../angle_harmonic.cpp
|
||||
rm ../angle_hybrid.cpp
|
||||
|
@ -100,12 +109,16 @@ elif (test $1 = 0) then
|
|||
rm ../improper_cvff.cpp
|
||||
rm ../improper_harmonic.cpp
|
||||
rm ../improper_hybrid.cpp
|
||||
rm ../improper_umbrella.cpp
|
||||
rm ../pair_hbond_dreiding_lj.cpp
|
||||
rm ../pair_hbond_dreiding_morse.cpp
|
||||
rm ../pair_lj_charmm_coul_charmm.cpp
|
||||
rm ../pair_lj_charmm_coul_charmm_implicit.cpp
|
||||
|
||||
rm ../angle_charmm.h
|
||||
rm ../angle_cosine.h
|
||||
rm ../angle_cosine_delta.h
|
||||
rm ../angle_cosine_periodic.h
|
||||
rm ../angle_cosine_squared.h
|
||||
rm ../angle_harmonic.h
|
||||
rm ../angle_hybrid.h
|
||||
|
@ -133,6 +146,9 @@ elif (test $1 = 0) then
|
|||
rm ../improper_cvff.h
|
||||
rm ../improper_harmonic.h
|
||||
rm ../improper_hybrid.h
|
||||
rm ../improper_umbrella.h
|
||||
rm ../pair_hbond_dreiding_lj.h
|
||||
rm ../pair_hbond_dreiding_morse.h
|
||||
rm ../pair_lj_charmm_coul_charmm.h
|
||||
rm ../pair_lj_charmm_coul_charmm_implicit.h
|
||||
|
||||
|
|
|
@ -0,0 +1,290 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Tod A Pascal (Caltech)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "angle_cosine_periodic.h"
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define SMALL 0.001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosinePeriodic::~AngleCosinePeriodic()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->sfree(setflag);
|
||||
memory->sfree(k);
|
||||
memory->sfree(b);
|
||||
memory->sfree(multiplicity);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosinePeriodic::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,i1,i2,i3,n,m,type,b_factor;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
|
||||
double tn,tn_1,tn_2,un,un_1,un_2;
|
||||
|
||||
eangle = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int **anglelist = neighbor->anglelist;
|
||||
int nanglelist = neighbor->nanglelist;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
for (n = 0; n < nanglelist; n++) {
|
||||
i1 = anglelist[n][0];
|
||||
i2 = anglelist[n][1];
|
||||
i3 = anglelist[n][2];
|
||||
type = anglelist[n][3];
|
||||
|
||||
// 1st bond
|
||||
|
||||
delx1 = x[i1][0] - x[i2][0];
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
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];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// c = cosine of angle
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
m = multiplicity[type];
|
||||
b_factor = b[type];
|
||||
|
||||
// cos(n*x) = Tn(cos(x))
|
||||
// Tn(x) = Chebyshev polynomials of the first kind: T_0 = 1, T_1 = x, ...
|
||||
// recurrence relationship:
|
||||
// Tn(x) = 2*x*T[n-1](x) - T[n-2](x) where T[-1](x) = 0
|
||||
// also, dTn(x)/dx = n*U[n-1](x)
|
||||
// where Un(x) = 2*x*U[n-1](x) - U[n-2](x) and U[-1](x) = 0
|
||||
// finally need to handle special case for n = 1
|
||||
|
||||
tn = 1.0;
|
||||
tn_1 = 1.0;
|
||||
tn_2 = 0.0;
|
||||
un = 1.0;
|
||||
un_1 = 2.0;
|
||||
un_2 = 0.0;
|
||||
|
||||
s = sqrt(1.0 - c*c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
s = 1.0/s;
|
||||
|
||||
// force & energy
|
||||
|
||||
tn_2 = c;
|
||||
for (i = 1; i <= m; i++) {
|
||||
tn = 2*c*tn_1 - tn_2;
|
||||
tn_2 = tn_1;
|
||||
tn_1 = tn;
|
||||
}
|
||||
|
||||
for (i = 2; i <= m; i++) {
|
||||
un = 2*c*un_1 - un_2;
|
||||
un_2 = un_1;
|
||||
un_1 = un;
|
||||
}
|
||||
tn = b_factor*pow((-1),m)*tn;
|
||||
un = b_factor*pow((-1),m)*m*un;
|
||||
|
||||
if (eflag) eangle = 2*k[type]*(1.0 - tn);
|
||||
|
||||
a = -k[type]*un;
|
||||
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;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
if (newton_bond || i1 < nlocal) {
|
||||
f[i1][0] += f1[0];
|
||||
f[i1][1] += f1[1];
|
||||
f[i1][2] += f1[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i2 < nlocal) {
|
||||
f[i2][0] -= f1[0] + f3[0];
|
||||
f[i2][1] -= f1[1] + f3[1];
|
||||
f[i2][2] -= f1[2] + f3[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i3 < nlocal) {
|
||||
f[i3][0] += f3[0];
|
||||
f[i3][1] += f3[1];
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosinePeriodic::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
|
||||
k = (double *) memory->smalloc((n+1)*sizeof(double),"angle:k");
|
||||
multiplicity = (int *) memory->smalloc((n+1)*sizeof(int),
|
||||
"angle:multiplicity");
|
||||
b = (int *) memory->smalloc((n+1)*sizeof(int),"angle:b");
|
||||
|
||||
setflag = (int *) memory->smalloc((n+1)*sizeof(int),"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more types
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosinePeriodic::coeff(int which, int narg, char **arg)
|
||||
{
|
||||
if (which > 0) return;
|
||||
if (narg != 4) error->all("Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
force->bounds(arg[0],atom->nangletypes,ilo,ihi);
|
||||
|
||||
double c_one = atof(arg[1]);
|
||||
int b_one = atoi(arg[2]);
|
||||
int n_one = atoi(arg[3]);
|
||||
if (n_one <= 0) error->all("Incorrect args for angle coefficients");
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = c_one/(n_one*n_one);
|
||||
b[i] = b_one;
|
||||
multiplicity[i] = n_one;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all("Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleCosinePeriodic::equilibrium_angle(int i)
|
||||
{
|
||||
return PI;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosinePeriodic::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&b[1],sizeof(int),atom->nangletypes,fp);
|
||||
fwrite(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosinePeriodic::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&b[1],sizeof(int),atom->nangletypes,fp);
|
||||
fread(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&b[1],atom->nangletypes,MPI_INT,0,world);
|
||||
MPI_Bcast(&multiplicity[1],atom->nangletypes,MPI_INT,0,world);
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleCosinePeriodic::single(int type, int i1, int i2, int i3)
|
||||
{
|
||||
double **x = atom->x;
|
||||
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
c = cos(acos(c)*multiplicity[type]);
|
||||
return k[type]*(1.0-b[type]*pow(-1,multiplicity[type])*c);
|
||||
}
|
|
@ -0,0 +1,49 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(cosine/periodic, AngleCosinePeriodic)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_ANGLE_PERIODIC_H
|
||||
#define LMP_ANGLE_PERIODIC_H
|
||||
|
||||
#include "stdio.h"
|
||||
#include "angle.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class AngleCosinePeriodic : public Angle {
|
||||
public:
|
||||
AngleCosinePeriodic(class LAMMPS *);
|
||||
~AngleCosinePeriodic();
|
||||
void compute(int, int);
|
||||
void coeff(int, int, char **);
|
||||
double equilibrium_angle(int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
double single(int, int, int, int);
|
||||
|
||||
private:
|
||||
double *k;
|
||||
int *multiplicity,*b;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
|
@ -0,0 +1,309 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Tod A Pascal (Caltech)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "mpi.h"
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "improper_umbrella.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define TOLERANCE 0.05
|
||||
#define SMALL 0.001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ImproperUmbrella::ImproperUmbrella(LAMMPS *lmp) : Improper(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ImproperUmbrella::~ImproperUmbrella()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->sfree(setflag);
|
||||
memory->sfree(kw);
|
||||
memory->sfree(w0);
|
||||
memory->sfree(C);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ImproperUmbrella::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,i4,n,type;
|
||||
double eimproper,f1[3],f2[3],f3[3],f4[3];
|
||||
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
|
||||
double domega,c,a,s,projhfg,dhax,dhay,dhaz,dahx,dahy,dahz,cotphi;
|
||||
double ax,ay,az,ra2,rh2,ra,rh,rar,rhr,arx,ary,arz,hrx,hry,hrz;
|
||||
|
||||
eimproper = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int **improperlist = neighbor->improperlist;
|
||||
int nimproperlist = neighbor->nimproperlist;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
for (n = 0; n < nimproperlist; n++) {
|
||||
i1 = improperlist[n][0];
|
||||
i2 = improperlist[n][1];
|
||||
i3 = improperlist[n][2];
|
||||
i4 = improperlist[n][3];
|
||||
type = improperlist[n][4];
|
||||
|
||||
// 1st bond
|
||||
|
||||
vb1x = x[i2][0] - x[i1][0];
|
||||
vb1y = x[i2][1] - x[i1][1];
|
||||
vb1z = x[i2][2] - x[i1][2];
|
||||
domain->minimum_image(vb1x,vb1y,vb1z);
|
||||
|
||||
// 2nd bond
|
||||
|
||||
vb2x = x[i3][0] - x[i1][0];
|
||||
vb2y = x[i3][1] - x[i1][1];
|
||||
vb2z = x[i3][2] - x[i1][2];
|
||||
domain->minimum_image(vb2x,vb2y,vb2z);
|
||||
|
||||
// 3rd bond
|
||||
|
||||
vb3x = x[i4][0] - x[i1][0];
|
||||
vb3y = x[i4][1] - x[i1][1];
|
||||
vb3z = x[i4][2] - x[i1][2];
|
||||
domain->minimum_image(vb3x,vb3y,vb3z);
|
||||
|
||||
// c0 calculation
|
||||
// A = vb1 X vb2 is perpendicular to IJK plane
|
||||
|
||||
ax = vb1y*vb2z-vb1z*vb2y;
|
||||
ay = vb1z*vb2x-vb1x*vb2z;
|
||||
az = vb1x*vb2y-vb1y*vb2x;
|
||||
ra2 = ax*ax+ay*ay+az*az;
|
||||
rh2 = vb3x*vb3x+vb3y*vb3y+vb3z*vb3z;
|
||||
ra = sqrt(ra2);
|
||||
rh = sqrt(rh2);
|
||||
if (ra < SMALL) ra = SMALL;
|
||||
if (rh < SMALL) rh = SMALL;
|
||||
|
||||
rar = 1/ra;
|
||||
rhr = 1/rh;
|
||||
arx = ax*rar;
|
||||
ary = ay*rar;
|
||||
arz = az*rar;
|
||||
hrx = vb3x*rhr;
|
||||
hry = vb3y*rhr;
|
||||
hrz = vb3z*rhr;
|
||||
|
||||
c = arx*hrx+ary*hry+arz*hrz;
|
||||
|
||||
// error check
|
||||
|
||||
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
|
||||
int me;
|
||||
MPI_Comm_rank(world,&me);
|
||||
if (screen) {
|
||||
fprintf(screen,"Improper problem: %d %d %d %d %d %d\n",
|
||||
me,update->ntimestep,
|
||||
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
|
||||
fprintf(screen," 1st atom: %d %g %g %g\n",
|
||||
me,x[i1][0],x[i1][1],x[i1][2]);
|
||||
fprintf(screen," 2nd atom: %d %g %g %g\n",
|
||||
me,x[i2][0],x[i2][1],x[i2][2]);
|
||||
fprintf(screen," 3rd atom: %d %g %g %g\n",
|
||||
me,x[i3][0],x[i3][1],x[i3][2]);
|
||||
fprintf(screen," 4th atom: %d %g %g %g\n",
|
||||
me,x[i4][0],x[i4][1],x[i4][2]);
|
||||
}
|
||||
}
|
||||
|
||||
if (c > 1.0) s = 1.0;
|
||||
if (c < -1.0) s = -1.0;
|
||||
|
||||
s = sqrt(1.0 - c*c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
cotphi = c/s;
|
||||
|
||||
projhfg = (vb3x*vb1x+vb3y*vb1y+vb3z*vb1z) /
|
||||
sqrt(vb1x*vb1x+vb1y*vb1y+vb1z*vb1z);
|
||||
projhfg += (vb3x*vb2x+vb3y*vb2y+vb3z*vb2z) /
|
||||
sqrt(vb2x*vb2x+vb2y*vb2y+vb2z*vb2z);
|
||||
if (projhfg > 0.0) {
|
||||
s *= -1.0;
|
||||
cotphi *= -1.0;
|
||||
}
|
||||
|
||||
// force and energy
|
||||
// if w0 = 0: E = k * (1 - cos w)
|
||||
// if w0 != 0: E = 0.5 * C (cos w - cos w0)^2, C = k/(sin(w0)^2
|
||||
|
||||
if (w0[type] == 0.0) {
|
||||
if (eflag) eimproper = kw[type] * (1.0-s);
|
||||
a = -kw[type];
|
||||
} else {
|
||||
domega = s - cos(w0[type]);
|
||||
a = 0.5 * C[type] * domega;
|
||||
if (eflag) eimproper = a * domega;
|
||||
a *= 2.0;
|
||||
}
|
||||
|
||||
// dhax = diffrence between H and A in X direction, etc
|
||||
|
||||
a = a*cotphi;
|
||||
dhax = hrx-c*arx;
|
||||
dhay = hry-c*ary;
|
||||
dhaz = hrz-c*arz;
|
||||
|
||||
dahx = arx-c*hrx;
|
||||
dahy = ary-c*hry;
|
||||
dahz = arz-c*hrz;
|
||||
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
|
||||
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
|
||||
|
||||
f4[0] = dahx*rhr;
|
||||
f4[1] = dahy*rhr;
|
||||
f4[2] = dahz*rhr;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
f1[2] = -(f2[2] + f3[2] + f4[2]);
|
||||
|
||||
// apply force to each of 4 atoms
|
||||
|
||||
if (newton_bond || i1 < nlocal) {
|
||||
f[i1][0] += f1[0]*a;
|
||||
f[i1][1] += f1[1]*a;
|
||||
f[i1][2] += f1[2]*a;
|
||||
}
|
||||
|
||||
if (newton_bond || i2 < nlocal) {
|
||||
f[i2][0] += f3[0]*a;
|
||||
f[i2][1] += f3[1]*a;
|
||||
f[i2][2] += f3[2]*a;
|
||||
}
|
||||
|
||||
if (newton_bond || i3 < nlocal) {
|
||||
f[i3][0] += f2[0]*a;
|
||||
f[i3][1] += f2[1]*a;
|
||||
f[i3][2] += f2[2]*a;
|
||||
}
|
||||
|
||||
if (newton_bond || i4 < nlocal) {
|
||||
f[i4][0] += f4[0]*a;
|
||||
f[i4][1] += f4[1]*a;
|
||||
f[i4][2] += f4[2]*a;
|
||||
}
|
||||
|
||||
if (evflag)
|
||||
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
|
||||
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ImproperUmbrella::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nimpropertypes;
|
||||
|
||||
kw = (double *) memory->smalloc((n+1)*sizeof(double),"improper:kw");
|
||||
w0 = (double *) memory->smalloc((n+1)*sizeof(double),"improper:w0");
|
||||
C = (double *) memory->smalloc((n+1)*sizeof(double),"improper:C");
|
||||
|
||||
setflag = (int *) memory->smalloc((n+1)*sizeof(int),"improper:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one type
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperUmbrella::coeff(int which, int narg, char **arg)
|
||||
{
|
||||
if (which > 0) return;
|
||||
if (narg != 3) error->all("Incorrect args for improper coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
force->bounds(arg[0],atom->nimpropertypes,ilo,ihi);
|
||||
|
||||
double k_one = atof(arg[1]);
|
||||
double w_one = atof(arg[2]);
|
||||
|
||||
// convert w0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
kw[i] = k_one;
|
||||
w0[i] = w_one/180.0 * PI;
|
||||
if (w_one == 0) C[i] = 1.0;
|
||||
else C[i] = kw[i]/(pow(sin(w0[i]),2));
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all("Incorrect args for improper coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperUmbrella::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&kw[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fwrite(&w0[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fwrite(&C[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperUmbrella::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&kw[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fread(&w0[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fread(&C[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
}
|
||||
MPI_Bcast(&kw[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&w0[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&C[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
|
||||
|
||||
for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
|
||||
}
|
|
@ -0,0 +1,46 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 IMPROPER_CLASS
|
||||
|
||||
ImproperStyle(umbrella,ImproperUmbrella)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_IMPROPER_UMBRELLA_H
|
||||
#define LMP_IMPROPER_UMBRELLA_H
|
||||
|
||||
#include "stdio.h"
|
||||
#include "improper.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ImproperUmbrella : public Improper {
|
||||
public:
|
||||
ImproperUmbrella(class LAMMPS *);
|
||||
~ImproperUmbrella();
|
||||
void compute(int, int);
|
||||
void coeff(int, int, char **);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
|
||||
private:
|
||||
double *kw, *w0, *C;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
|
@ -0,0 +1,417 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Tod A Pascal (Caltech)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_hbond_dreiding_lj.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "domain.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
||||
#define MAX(a,b) ((a) > (b) ? (a) : (b))
|
||||
|
||||
#define SMALL 0.001
|
||||
#define CHUNK 8
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
single_enable = 0;
|
||||
|
||||
// hbond cannot compute virial as F dot r
|
||||
// due to using map() to find bonded H atoms which are not near donor atom
|
||||
|
||||
no_virial_compute = 1;
|
||||
|
||||
PI = 4.0*atan(1.0);
|
||||
|
||||
nparams = maxparam = 0;
|
||||
params = NULL;
|
||||
|
||||
nextra = 1;
|
||||
pvector = new double[1];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairHbondDreidingLJ::~PairHbondDreidingLJ()
|
||||
{
|
||||
memory->sfree(params);
|
||||
delete [] pvector;
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy_2d_int_array(setflag);
|
||||
memory->destroy_2d_double_array(cutsq);
|
||||
|
||||
delete [] donor;
|
||||
delete [] acceptor;
|
||||
memory->destroy_3d_int_array(type2param);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingLJ::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
|
||||
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
|
||||
double factor_hb,force_angle,force_kernel,evdwl;
|
||||
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
|
||||
double fi[3],fj[3],delr1[3],delr2[3];
|
||||
double r2inv,r10inv;
|
||||
int *ilist,*jlist,*klist,*numneigh,**firstneigh;
|
||||
Param *pm;
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int **special = atom->special;
|
||||
int *type = atom->type;
|
||||
int **nspecial = atom->nspecial;
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
double *special_lj = force->special_lj;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// ii = loop over donors
|
||||
// jj = loop over acceptors
|
||||
// kk = loop over hydrogens bonded to donor
|
||||
|
||||
int hbcount = 0;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itype = type[i];
|
||||
if (!donor[itype]) continue;
|
||||
klist = special[i];
|
||||
knum = nspecial[i][0];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
if (j < nall) factor_hb = 1.0;
|
||||
else {
|
||||
factor_hb = special_lj[j/nall];
|
||||
j %= nall;
|
||||
}
|
||||
jtype = type[j];
|
||||
if (!acceptor[jtype]) continue;
|
||||
|
||||
delx = x[i][0] - x[j][0];
|
||||
dely = x[i][1] - x[j][1];
|
||||
delz = x[i][2] - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
|
||||
for (kk = 0; kk < knum; kk++) {
|
||||
k = atom->map(klist[kk]);
|
||||
if (k < 0) continue;
|
||||
ktype = type[k];
|
||||
m = type2param[itype][jtype][ktype];
|
||||
if (m < 0) continue;
|
||||
pm = ¶ms[m];
|
||||
|
||||
if (rsq < pm->cutsq) {
|
||||
delr1[0] = x[i][0] - x[k][0];
|
||||
delr1[1] = x[i][1] - x[k][1];
|
||||
delr1[2] = x[i][2] - x[k][2];
|
||||
domain->minimum_image(delr1);
|
||||
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
delr2[0] = x[j][0] - x[k][0];
|
||||
delr2[1] = x[j][1] - x[k][1];
|
||||
delr2[2] = x[j][2] - x[k][2];
|
||||
domain->minimum_image(delr2);
|
||||
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2];
|
||||
c /= r1*r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
ac = acos(c);
|
||||
|
||||
if (ac > pm->cut_angle && ac < (2.0*PI - pm->cut_angle)) {
|
||||
s = sqrt(1.0 - c*c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
|
||||
// LJ-specific kernel
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r10inv = r2inv*r2inv*r2inv*r2inv*r2inv;
|
||||
force_kernel = r10inv*(pm->lj1*r2inv - pm->lj2)*r2inv *
|
||||
pow(c,pm->ap);
|
||||
force_angle = pm->ap * r10inv*(pm->lj3*r2inv - pm->lj4) *
|
||||
pow(c,pm->ap-1)*s;
|
||||
|
||||
if (eflag) {
|
||||
evdwl = r10inv*(pm->lj3*r2inv - pm->lj4) * pow(c,pm->ap) -
|
||||
pm->offset;
|
||||
evdwl *= factor_hb;
|
||||
}
|
||||
|
||||
a = factor_hb*force_angle/s;
|
||||
b = factor_hb*force_kernel;
|
||||
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
|
||||
vx1 = a11*delr1[0] + a12*delr2[0];
|
||||
vx2 = a22*delr2[0] + a12*delr1[0];
|
||||
vy1 = a11*delr1[1] + a12*delr2[1];
|
||||
vy2 = a22*delr2[1] + a12*delr1[1];
|
||||
vz1 = a11*delr1[2] + a12*delr2[2];
|
||||
vz2 = a22*delr2[2] + a12*delr1[2];
|
||||
|
||||
fi[0] = vx1 + b*delx;
|
||||
fi[1] = vy1 + b*dely;
|
||||
fi[2] = vz1 + b*delz;
|
||||
fj[0] = vx2 - b*delx;
|
||||
fj[1] = vy2 - b*dely;
|
||||
fj[2] = vz2 - b*delz;
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
|
||||
f[k][0] -= vx1 + vx2;
|
||||
f[k][1] -= vy1 + vy2;
|
||||
f[k][2] -= vz1 + vz2;
|
||||
|
||||
// KIJ instead of IJK b/c delr1/delr2 are both with respect to k
|
||||
|
||||
if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2);
|
||||
|
||||
hbcount++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (eflag_global) pvector[0] = hbcount;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingLJ::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
// mark all setflag as set, since don't require pair_coeff of all I,J
|
||||
|
||||
setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 1;
|
||||
|
||||
cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
|
||||
|
||||
donor = new int[n+1];
|
||||
acceptor = new int[n+1];
|
||||
type2param = memory->create_3d_int_array(n+1,n+1,n+1,"pair:type2param");
|
||||
|
||||
int i,j,k;
|
||||
for (i = 1; i <= n; i++)
|
||||
for (j = 1; j <= n; j++)
|
||||
for (k = 1; k <= n; k++)
|
||||
type2param[i][j][k] = -1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingLJ::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all("Illegal pair_style command");
|
||||
|
||||
ap_global = force->inumeric(arg[0]);
|
||||
cut_global = force->numeric(arg[1]);
|
||||
cut_angle_global = force->numeric(arg[2]) * PI/180.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingLJ::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 6 || narg > 9)
|
||||
error->all("Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi,klo,khi;
|
||||
force->bounds(arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(arg[1],atom->ntypes,jlo,jhi);
|
||||
force->bounds(arg[2],atom->ntypes,klo,khi);
|
||||
|
||||
int donor_flag;
|
||||
if (strcmp(arg[3],"i") == 0) donor_flag = 0;
|
||||
else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
|
||||
else error->all("Incorrect args for pair coefficients");
|
||||
|
||||
double epsilon_one = force->numeric(arg[4]);
|
||||
double sigma_one = force->numeric(arg[5]);
|
||||
|
||||
int ap_one = ap_global;
|
||||
if (narg == 7) ap_one = force->inumeric(arg[6]);
|
||||
double cut_one = cut_global;
|
||||
if (narg == 8) cut_one = force->numeric(arg[7]);
|
||||
double cut_angle_one = cut_angle_global;
|
||||
if (narg == 9) cut_angle_one = force->numeric(arg[8]) * PI/180.0;
|
||||
|
||||
// grow params array if necessary
|
||||
|
||||
if (nparams == maxparam) {
|
||||
maxparam += CHUNK;
|
||||
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
|
||||
"pair:params");
|
||||
}
|
||||
|
||||
params[nparams].epsilon = epsilon_one;
|
||||
params[nparams].sigma = sigma_one;
|
||||
params[nparams].ap = ap_one;
|
||||
params[nparams].cut = cut_one;
|
||||
params[nparams].cutsq = cut_one*cut_one;
|
||||
params[nparams].cut_angle = cut_angle_one;
|
||||
|
||||
// flag type2param with either i,j = D,A or j,i = D,A
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++)
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++)
|
||||
for (int k = klo; k <= khi; k++) {
|
||||
if (donor_flag == 0) type2param[i][j][k] = nparams;
|
||||
else type2param[j][i][k] = nparams;
|
||||
count++;
|
||||
}
|
||||
nparams++;
|
||||
|
||||
if (count == 0) error->all("Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingLJ::init_style()
|
||||
{
|
||||
// molecular system required to use special list to find H atoms
|
||||
// tags required to use special list
|
||||
// pair newton on required since are looping over D atoms
|
||||
// and computing forces on A,H which may be on different procs
|
||||
|
||||
if (atom->molecular == 0)
|
||||
error->all("Pair style hbond/dreiding requires molecular system");
|
||||
if (atom->tag_enable == 0)
|
||||
error->all("Pair style hbond/dreiding requires atom IDs");
|
||||
if (atom->map_style == 0)
|
||||
error->all("Pair style hbond/dreiding requires an atom map, "
|
||||
"see atom_modify");
|
||||
if (force->newton_pair == 0)
|
||||
error->all("Pair style hbond/dreiding requires newton pair on");
|
||||
|
||||
// set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor
|
||||
|
||||
int anyflag = 0;
|
||||
int n = atom->ntypes;
|
||||
for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0;
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = 1; j <= n; j++)
|
||||
for (int k = 1; k <= n; k++)
|
||||
if (type2param[i][j][k] >= 0) {
|
||||
anyflag = 1;
|
||||
donor[i] = 1;
|
||||
acceptor[j] = 1;
|
||||
}
|
||||
|
||||
if (!anyflag) error->all("No pair hbond/dreiding coefficients set");
|
||||
|
||||
// set additional param values
|
||||
// offset is for LJ only, angle term is not included
|
||||
|
||||
for (int m = 0; m < nparams; m++) {
|
||||
params[m].lj1 = 60.0*params[m].epsilon*pow(params[m].sigma,12.0);
|
||||
params[m].lj2 = 60.0*params[m].epsilon*pow(params[m].sigma,10.0);
|
||||
params[m].lj3 = 5.0*params[m].epsilon*pow(params[m].sigma,12.0);
|
||||
params[m].lj4 = 6.0*params[m].epsilon*pow(params[m].sigma,10.0);
|
||||
|
||||
if (offset_flag) {
|
||||
double ratio = params[m].sigma / params[m].cut;
|
||||
params[m].offset = params[m].epsilon *
|
||||
((2.0*pow(ratio,9.0)) - (3.0*pow(ratio,6.0)));
|
||||
} else params[m].offset = 0.0;
|
||||
}
|
||||
|
||||
// full neighbor list request
|
||||
|
||||
int irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairHbondDreidingLJ::init_one(int i, int j)
|
||||
{
|
||||
int m;
|
||||
|
||||
// return maximum cutoff for any K with I,J = D,A or J,I = D,A
|
||||
// donor/acceptor is not symmetric, IJ interaction != JI interaction
|
||||
|
||||
double cut = 0.0;
|
||||
for (int k = 1; k <= atom->ntypes; k++) {
|
||||
m = type2param[i][j][k];
|
||||
if (m >= 0) cut = MAX(cut,params[m].cut);
|
||||
m = type2param[j][i][k];
|
||||
if (m >= 0) cut = MAX(cut,params[m].cut);
|
||||
}
|
||||
return cut;
|
||||
}
|
|
@ -0,0 +1,65 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_HBOND_DREIDING_LJ_H
|
||||
#define LMP_PAIR_HBOND_DREIDING_LJ_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairHbondDreidingLJ : public Pair {
|
||||
public:
|
||||
PairHbondDreidingLJ(class LAMMPS *);
|
||||
virtual ~PairHbondDreidingLJ();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
virtual void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
protected:
|
||||
double cut_global,cut_angle_global;
|
||||
int ap_global;
|
||||
double PI;
|
||||
|
||||
struct Param {
|
||||
double epsilon,sigma;
|
||||
double lj1,lj2,lj3,lj4;
|
||||
double d0,alpha,r0;
|
||||
double morse1;
|
||||
double cut,cut_angle,cutsq,offset;
|
||||
int ap;
|
||||
};
|
||||
|
||||
Param *params; // parameter set for an I-J-K interaction
|
||||
int nparams; // number of parameters read
|
||||
int maxparam;
|
||||
|
||||
int *donor; // 1 if this type is ever a donor, else 0
|
||||
int *acceptor; // 1 if this type is ever an acceptor, else 0
|
||||
int ***type2param; // mapping from D,A,H to params, -1 if no map
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
|
@ -0,0 +1,318 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Tod A Pascal (Caltech)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_hbond_dreiding_morse.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "domain.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
||||
#define MAX(a,b) ((a) > (b) ? (a) : (b))
|
||||
|
||||
#define SMALL 0.001
|
||||
#define CHUNK 8
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairHbondDreidingMorse::PairHbondDreidingMorse(LAMMPS *lmp) :
|
||||
PairHbondDreidingLJ(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingMorse::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
|
||||
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
|
||||
double factor_hb,force_angle,force_kernel,evdwl;
|
||||
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
|
||||
double fi[3],fj[3],delr1[3],delr2[3];
|
||||
double r,dr,dexp,emorse;
|
||||
int *ilist,*jlist,*klist,*numneigh,**firstneigh;
|
||||
Param *pm;
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int **special = atom->special;
|
||||
int *type = atom->type;
|
||||
int **nspecial = atom->nspecial;
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
double *special_lj = force->special_lj;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// ii = loop over donors
|
||||
// jj = loop over acceptors
|
||||
// kk = loop over hydrogens bonded to donor
|
||||
|
||||
int hbcount = 0;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itype = type[i];
|
||||
if (!donor[itype]) continue;
|
||||
klist = special[i];
|
||||
knum = nspecial[i][0];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
if (j < nall) factor_hb = 1.0;
|
||||
else {
|
||||
factor_hb = special_lj[j/nall];
|
||||
j %= nall;
|
||||
}
|
||||
jtype = type[j];
|
||||
if (!acceptor[jtype]) continue;
|
||||
|
||||
delx = x[i][0] - x[j][0];
|
||||
dely = x[i][1] - x[j][1];
|
||||
delz = x[i][2] - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
|
||||
for (kk = 0; kk < knum; kk++) {
|
||||
k = atom->map(klist[kk]);
|
||||
if (k < 0) continue;
|
||||
ktype = type[k];
|
||||
m = type2param[itype][jtype][ktype];
|
||||
if (m < 0) continue;
|
||||
pm = ¶ms[m];
|
||||
|
||||
if (rsq < pm->cutsq) {
|
||||
delr1[0] = x[i][0] - x[k][0];
|
||||
delr1[1] = x[i][1] - x[k][1];
|
||||
delr1[2] = x[i][2] - x[k][2];
|
||||
domain->minimum_image(delr1);
|
||||
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
delr2[0] = x[j][0] - x[k][0];
|
||||
delr2[1] = x[j][1] - x[k][1];
|
||||
delr2[2] = x[j][2] - x[k][2];
|
||||
domain->minimum_image(delr2);
|
||||
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2];
|
||||
c /= r1*r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
ac = acos(c);
|
||||
|
||||
if (ac > pm->cut_angle && ac < (2.0*PI - pm->cut_angle)) {
|
||||
s = sqrt(1.0 - c*c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
|
||||
// Morse-specific kernel
|
||||
|
||||
r = sqrt(rsq);
|
||||
dr = r - pm->r0;
|
||||
dexp = exp(-pm->alpha * dr);
|
||||
force_kernel = pm->morse1*(dexp*dexp - dexp)/r * pow(c,pm->ap);
|
||||
emorse = pm->d0 * (dexp*dexp - 2.0*dexp) - pm->offset;
|
||||
force_angle = pm->ap * emorse * pow(c,pm->ap-1)*s;
|
||||
|
||||
if (eflag) {
|
||||
evdwl = emorse * pow(c,params[m].ap);
|
||||
evdwl *= factor_hb;
|
||||
}
|
||||
|
||||
a = factor_hb*force_angle/s;
|
||||
b = factor_hb*force_kernel;
|
||||
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
|
||||
vx1 = a11*delr1[0] + a12*delr2[0];
|
||||
vx2 = a22*delr2[0] + a12*delr1[0];
|
||||
vy1 = a11*delr1[1] + a12*delr2[1];
|
||||
vy2 = a22*delr2[1] + a12*delr1[1];
|
||||
vz1 = a11*delr1[2] + a12*delr2[2];
|
||||
vz2 = a22*delr2[2] + a12*delr1[2];
|
||||
|
||||
fi[0] = vx1 + b*delx;
|
||||
fi[1] = vy1 + b*dely;
|
||||
fi[2] = vz1 + b*delz;
|
||||
fj[0] = vx2 - b*delx;
|
||||
fj[1] = vy2 - b*dely;
|
||||
fj[2] = vz2 - b*delz;
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
|
||||
f[k][0] -= vx1 + vx2;
|
||||
f[k][1] -= vy1 + vy2;
|
||||
f[k][2] -= vz1 + vz2;
|
||||
|
||||
// KIJ instead of IJK b/c delr1/delr2 are both with respect to k
|
||||
|
||||
if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2);
|
||||
|
||||
hbcount++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingMorse::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 7 || narg > 10)
|
||||
error->all("Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi,klo,khi;
|
||||
force->bounds(arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(arg[1],atom->ntypes,jlo,jhi);
|
||||
force->bounds(arg[2],atom->ntypes,klo,khi);
|
||||
|
||||
int donor_flag;
|
||||
if (strcmp(arg[3],"i") == 0) donor_flag = 0;
|
||||
else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
|
||||
else error->all("Incorrect args for pair coefficients");
|
||||
|
||||
double d0_one = force->numeric(arg[4]);
|
||||
double alpha_one = force->numeric(arg[5]);
|
||||
double r0_one = force->numeric(arg[6]);
|
||||
|
||||
int ap_one = ap_global;
|
||||
if (narg == 8) ap_one = force->inumeric(arg[7]);
|
||||
double cut_one = cut_global;
|
||||
if (narg == 9) cut_one = force->numeric(arg[8]);
|
||||
double cut_angle_one = cut_angle_global;
|
||||
if (narg == 10) cut_angle_one = force->numeric(arg[9]) * PI/180.0;
|
||||
|
||||
// grow params array if necessary
|
||||
|
||||
if (nparams == maxparam) {
|
||||
maxparam += CHUNK;
|
||||
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
|
||||
"pair:params");
|
||||
}
|
||||
|
||||
params[nparams].d0 = d0_one;
|
||||
params[nparams].alpha = alpha_one;
|
||||
params[nparams].r0 = r0_one;
|
||||
params[nparams].ap = ap_one;
|
||||
params[nparams].cut = cut_one;
|
||||
params[nparams].cutsq = cut_one*cut_one;
|
||||
params[nparams].cut_angle = cut_angle_one;
|
||||
|
||||
// flag type2param with either i,j = D,A or j,i = D,A
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++)
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++)
|
||||
for (int k = klo; k <= khi; k++) {
|
||||
if (donor_flag == 0) type2param[i][j][k] = nparams;
|
||||
else type2param[j][i][k] = nparams;
|
||||
count++;
|
||||
}
|
||||
nparams++;
|
||||
|
||||
if (count == 0) error->all("Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairHbondDreidingMorse::init_style()
|
||||
{
|
||||
// molecular system required to use special list to find H atoms
|
||||
// tags required to use special list
|
||||
// pair newton on required since are looping over D atoms
|
||||
// and computing forces on A,H which may be on different procs
|
||||
|
||||
if (atom->molecular == 0)
|
||||
error->all("Pair style hbond/dreiding requires molecular system");
|
||||
if (atom->tag_enable == 0)
|
||||
error->all("Pair style hbond/dreiding requires atom IDs");
|
||||
if (atom->map_style == 0)
|
||||
error->all("Pair style hbond/dreiding requires an atom map, "
|
||||
"see atom_modify");
|
||||
if (force->newton_pair == 0)
|
||||
error->all("Pair style hbond/dreiding requires newton pair on");
|
||||
|
||||
// set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor
|
||||
|
||||
int anyflag = 0;
|
||||
int n = atom->ntypes;
|
||||
for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0;
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = 1; j <= n; j++)
|
||||
for (int k = 1; k <= n; k++)
|
||||
if (type2param[i][j][k] >= 0) {
|
||||
anyflag = 1;
|
||||
donor[i] = 1;
|
||||
acceptor[j] = 1;
|
||||
}
|
||||
|
||||
if (!anyflag) error->all("No pair hbond/dreiding coefficients set");
|
||||
|
||||
// set additional param values
|
||||
// offset is for Morse only, angle term is not included
|
||||
|
||||
for (int m = 0; m < nparams; m++) {
|
||||
params[m].morse1 = 2.0*params[m].d0*params[m].alpha;
|
||||
|
||||
if (offset_flag) {
|
||||
double alpha_dr = -params[m].alpha * (params[m].cut - params[m].r0);
|
||||
params[m].offset = params[m].d0 *
|
||||
((exp(2.0*alpha_dr)) - (2.0*exp(alpha_dr)));
|
||||
} else params[m].offset = 0.0;
|
||||
}
|
||||
|
||||
// full neighbor list request
|
||||
|
||||
int irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
|
@ -0,0 +1,39 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(hbond/dreiding/morse,PairHbondDreidingMorse)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_H
|
||||
#define LMP_PAIR_HBOND_DREIDING_MORSE_H
|
||||
|
||||
#include "pair_hbond_dreiding_lj.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairHbondDreidingMorse : public PairHbondDreidingLJ {
|
||||
public:
|
||||
PairHbondDreidingMorse(class LAMMPS *);
|
||||
~PairHbondDreidingMorse() {}
|
||||
void compute(int, int);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
Loading…
Reference in New Issue