jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12434 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-09-08 15:44:38 +00:00
parent bcb5c5a6e5
commit a0b94e9053
11 changed files with 2258 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -15,7 +15,7 @@ OBJ = $(SRC:.cpp=.o)
PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
kokkos kspace manybody mc meam misc molecule mpiio opt peri poems \
reax replica rigid shock snap srd voronoi xtc
qeq reax replica rigid shock snap srd voronoi xtc
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
user-cuda user-eff user-fep user-intel user-lb user-misc \

795
src/QEQ/fix_qeq.cpp Normal file
View File

@ -0,0 +1,795 @@
/* ----------------------------------------------------------------------
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: Ray Shan (Sandia)
Based on fix qeq/reax by H. Metin Aktulga
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_qeq.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "kspace.h"
#include "group.h"
#include "pair.h"
#include "respa.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 8) error->all(FLERR,"Illegal fix qeq command");
nevery = force->inumeric(FLERR,arg[3]);
cutoff = force->numeric(FLERR,arg[4]);
tolerance = force->numeric(FLERR,arg[5]);
maxiter = force->inumeric(FLERR,arg[6]);
alpha = 0.20;
swa = 0.0;
swb = cutoff;
shld = NULL;
n = n_cap = 0;
N = nmax = 0;
m_fill = m_cap = 0;
pack_flag = 0;
s = NULL;
t = NULL;
nprev = 5;
Hdia_inv = NULL;
b_s = NULL;
b_t = NULL;
b_prc = NULL;
b_prm = NULL;
// CG
p = NULL;
q = NULL;
r = NULL;
d = NULL;
// H matrix
H.firstnbr = NULL;
H.numnbrs = NULL;
H.jlist = NULL;
H.val = NULL;
// others
cutoff_sq = cutoff*cutoff;
chizj = NULL;
qf = NULL;
q1 = NULL;
q2 = NULL;
comm_forward = comm_reverse = 1;
// perform initial allocation of atom-based arrays
// register with Atom class
s_hist = t_hist = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
for( int i = 0; i < atom->nmax; i++ )
for (int j = 0; j < nprev; ++j )
s_hist[i][j] = t_hist[i][j] = 0;
read_file(arg[7]);
}
/* ---------------------------------------------------------------------- */
FixQEq::~FixQEq()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
memory->destroy(s_hist);
memory->destroy(t_hist);
deallocate_storage();
deallocate_matrix();
memory->destroy(shld);
memory->destroy(chi);
memory->destroy(eta);
memory->destroy(gamma);
memory->destroy(zeta);
memory->destroy(zcore);
}
/* ---------------------------------------------------------------------- */
int FixQEq::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
mask |= MIN_PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixQEq::read_file(char *arg)
{
int i,itype,ntypes;
double v1,v2,v3,v4,v5;
char units[20];
FILE *pf;
ntypes = atom->ntypes;
memory->create(chi,ntypes+1,"qeq:chi");
memory->create(eta,ntypes+1,"qeq:eta");
memory->create(gamma,ntypes+1,"qeq:gamma");
memory->create(zeta,ntypes+1,"qeq:zeta");
memory->create(zcore,ntypes+1,"qeq:zcore");
if (comm->me == 0) {
if ((pf = fopen(arg,"r")) == NULL)
error->one(FLERR,"Fix qeq parameter file could not be found");
fscanf(pf,"%s",units);
if ((strcmp(units,"real") != 0) && (strcmp(units,"metal") != 0))
error->one(FLERR,"Fix qeq param file must be in real or metal units");
for (i = 1; i <= ntypes && !feof(pf); i++) {
fscanf(pf,"%d %lg %lg %lg %lg %lg",&itype,&v1,&v2,&v3,&v4,&v5);
if (itype < 1 || itype > ntypes)
error->one(FLERR,"Fix qeq invalid atom type in param file");
chi[itype] = v1;
eta[itype] = v2;
gamma[itype] = v3;
zeta[itype] = v4;
zcore[itype] = v5;
// unit conversion
if (strcmp(units,"real") == 0) {
if (strcmp(update->unit_style,"metal") == 0) {
chi[itype] *= 0.0433634;
eta[itype] *= 0.0433634;
gamma[itype] *= 1.0;
zeta[itype] *= 1.0;
zcore[itype] *= 1.0;
} else if (strcmp(update->unit_style,"si") == 0) {
chi[itype] *= 6.95e-21;
eta[itype] *= 6.95e-21;
gamma[itype] *= 1.0e-10;
zeta[itype] *= 1.0e10;
zcore[itype] *= 1.6021765e-19;
} else if (strcmp(update->unit_style,"cgs") == 0) {
chi[itype] *= 6.95e-14;
eta[itype] *= 6.95e-14;
gamma[itype] *= 1.0e-8;
zeta[itype] *= 1.0e8;
zcore[itype] *= 4.8032044e-10;
} else if (strcmp(update->unit_style,"electron") == 0) {
chi[itype] *= 0.00159362;
eta[itype] *= 0.00159362;
gamma[itype] *= 1.8903;
zeta[itype] *= 0.5290;
zcore[itype] *= 1.0;
}
} else if (strcmp(units,"metal") == 0) {
if (strcmp(update->unit_style,"real") == 0) {
chi[itype] *= 23.0609;
eta[itype] *= 23.0609;
gamma[itype] *= 1.0;
zeta[itype] *= 1.0;
zcore[itype] *= 1.0;
} else if (strcmp(update->unit_style,"si") == 0) {
chi[itype] *= 1.6021e-19;
eta[itype] *= 1.6021e-19;
gamma[itype] *= 1.0e-10;
zeta[itype] *= 1.0e10;
zcore[itype] *= 1.6021765e-19;
} else if (strcmp(update->unit_style,"cgs") == 0) {
chi[itype] *= 1.6021e-12;
eta[itype] *= 1.6021e-12;
gamma[itype] *= 1.0e-8;
zeta[itype] *= 1.0e8;
zcore[itype] *= 4.8032044e-10;
} else if (strcmp(update->unit_style,"electron") == 0) {
chi[itype] *= 0.0367502;
eta[itype] *= 0.0367502;
gamma[itype] *= 1.8903;
zeta[itype] *= 0.5290;
zcore[itype] *= 1.0;
}
}
}
if (i <= ntypes) error->one(FLERR,"Invalid param file for fix qeq");
fclose(pf);
}
MPI_Bcast(&chi[1],ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(&eta[1],ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(&gamma[1],ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(&zeta[1],ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(&zcore[1],ntypes+1,MPI_DOUBLE,0,world);
}
/* ---------------------------------------------------------------------- */
void FixQEq::allocate_storage()
{
nmax = atom->nmax;
memory->create(s,nmax,"qeq:s");
memory->create(t,nmax,"qeq:t");
memory->create(Hdia_inv,nmax,"qeq:Hdia_inv");
memory->create(b_s,nmax,"qeq:b_s");
memory->create(b_t,nmax,"qeq:b_t");
memory->create(b_prc,nmax,"qeq:b_prc");
memory->create(b_prm,nmax,"qeq:b_prm");
memory->create(p,nmax,"qeq:p");
memory->create(q,nmax,"qeq:q");
memory->create(r,nmax,"qeq:r");
memory->create(d,nmax,"qeq:d");
memory->create(chizj,nmax,"qeq:chizj");
memory->create(qf,nmax,"qeq:qf");
memory->create(q1,nmax,"qeq:q1");
memory->create(q2,nmax,"qeq:q2");
}
/* ---------------------------------------------------------------------- */
void FixQEq::deallocate_storage()
{
memory->destroy(s);
memory->destroy(t);
memory->destroy( Hdia_inv );
memory->destroy( b_s );
memory->destroy( b_t );
memory->destroy( b_prc );
memory->destroy( b_prm );
memory->destroy( p );
memory->destroy( q );
memory->destroy( r );
memory->destroy( d );
memory->destroy( chizj );
memory->destroy( qf );
memory->destroy( q1 );
memory->destroy( q2 );
}
/* ---------------------------------------------------------------------- */
void FixQEq::reallocate_storage()
{
deallocate_storage();
allocate_storage();
init_storage();
}
/* ---------------------------------------------------------------------- */
void FixQEq::allocate_matrix()
{
int i,ii,inum,m;
int *ilist, *numneigh;
int mincap;
double safezone;
mincap = MIN_CAP;
safezone = SAFE_ZONE;
n = atom->nlocal;
n_cap = MAX( (int)(n * safezone), mincap );
N = atom->nlocal + atom->nghost;
// determine the total space for the H matrix
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
m = 0;
for( ii = 0; ii < inum; ii++ ) {
i = ilist[ii];
m += numneigh[i];
}
m_cap = MAX( (int)(m * safezone), mincap * MIN_NBRS );
H.n = n_cap;
H.m = m_cap;
memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr");
memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs");
memory->create(H.jlist,m_cap,"qeq:H.jlist");
memory->create(H.val,m_cap,"qeq:H.val");
}
/* ---------------------------------------------------------------------- */
void FixQEq::deallocate_matrix()
{
memory->destroy( H.firstnbr );
memory->destroy( H.numnbrs );
memory->destroy( H.jlist );
memory->destroy( H.val );
}
/* ---------------------------------------------------------------------- */
void FixQEq::reallocate_matrix()
{
deallocate_matrix();
allocate_matrix();
}
/* ---------------------------------------------------------------------- */
void FixQEq::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void FixQEq::setup_pre_force(int vflag)
{
neighbor->build_one(list->index);
deallocate_storage();
allocate_storage();
init_storage();
deallocate_matrix();
allocate_matrix();
pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEq::setup_pre_force_respa(int vflag, int ilevel)
{
if (ilevel < nlevels_respa-1) return;
setup_pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEq::min_setup_pre_force(int vflag)
{
setup_pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEq::init_storage()
{
int N = atom->nlocal + atom->nghost;
for( int i = 0; i < N; i++ ) {
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
b_prc[i] = 0;
b_prm[i] = 0;
s[i] = t[i] = 0;
chizj[i] = 0.0;
qf[i] = 0.0;
q1[i] = 0.0;
q2[i] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
void FixQEq::pre_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixQEq::min_pre_force(int vflag)
{
pre_force(vflag);
}
/* ---------------------------------------------------------------------- */
int FixQEq::CG( double *b, double *x )
{
int i, j;
double tmp, alfa, beta, b_norm;
double sig_old, sig_new;
int nn, jj;
int *ilist;
nn = list->inum;
ilist = list->ilist;
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this ); //Coll_Vector( q );
vector_sum( r , 1., b, -1., q, nn );
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
d[j] = r[j] * Hdia_inv[j]; //pre-condition
}
b_norm = parallel_norm( b, nn );
sig_new = parallel_dot( r, d, nn);
for( i = 1; i < maxiter && sqrt(sig_new) / b_norm > tolerance; ++i ) {
comm->forward_comm_fix(this); //Dist_vector( d );
sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this); //Coll_vector( q );
tmp = parallel_dot( d, q, nn);
alfa = sig_new / tmp;
vector_add( x, alfa, d, nn );
vector_add( r, -alfa, q, nn );
// pre-conditioning
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
p[j] = r[j] * Hdia_inv[j];
}
sig_old = sig_new;
sig_new = parallel_dot( r, p, nn);
beta = sig_new / sig_old;
vector_sum( d, 1., p, beta, d, nn );
}
if (i >= maxiter && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq CG convergence failed (%g) after %d iterations "
"at " BIGINT_FORMAT " step",sqrt(sig_new) / b_norm,i,update->ntimestep);
error->warning(FLERR,str);
}
return i;
}
/* ---------------------------------------------------------------------- */
void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
nn = atom->nlocal;
NN = atom->nlocal + atom->nghost;
ilist = list->ilist;
for( i = 0; i < nn; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = eta[ atom->type[i] ] * x[i];
}
for( i = nn; i < NN; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = 0;
}
for( i = 0; i < nn; ++i ) {
if (atom->mask[i] & groupbit) {
for( itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixQEq::calculate_Q()
{
int i, k;
double u, s_sum, t_sum;
double *q = atom->q;
int nn, ii;
int *ilist;
nn = list->inum;
ilist = list->ilist;
s_sum = parallel_vector_acc( s, nn );
t_sum = parallel_vector_acc( t, nn);
u = s_sum / t_sum;
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
q[i] = s[i] - u * t[i];
/* backup s & t */
for( k = 4; k > 0; --k ) {
s_hist[i][k] = s_hist[i][k-1];
t_hist[i][k] = t_hist[i][k-1];
}
s_hist[i][0] = s[i];
t_hist[i][0] = t[i];
}
}
pack_flag = 4;
comm->forward_comm_fix( this ); //Dist_vector( atom->q );
}
/* ---------------------------------------------------------------------- */
int FixQEq::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int m;
if( pack_flag == 1)
for(m = 0; m < n; m++) buf[m] = d[list[m]];
else if( pack_flag == 2 )
for(m = 0; m < n; m++) buf[m] = s[list[m]];
else if( pack_flag == 3 )
for(m = 0; m < n; m++) buf[m] = t[list[m]];
else if( pack_flag == 4 )
for(m = 0; m < n; m++) buf[m] = atom->q[list[m]];
return n;
}
/* ---------------------------------------------------------------------- */
void FixQEq::unpack_forward_comm(int n, int first, double *buf)
{
int i, m;
if( pack_flag == 1)
for(m = 0, i = first; m < n; m++, i++) d[i] = buf[m];
else if( pack_flag == 2)
for(m = 0, i = first; m < n; m++, i++) s[i] = buf[m];
else if( pack_flag == 3)
for(m = 0, i = first; m < n; m++, i++) t[i] = buf[m];
else if( pack_flag == 4)
for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
}
/* ---------------------------------------------------------------------- */
int FixQEq::pack_reverse_comm(int n, int first, double *buf)
{
int i, m;
for(m = 0, i = first; m < n; m++, i++) buf[m] = q[i];
return n;
}
/* ---------------------------------------------------------------------- */
void FixQEq::unpack_reverse_comm(int n, int *list, double *buf)
{
int m;
for(m = 0; m < n; m++) q[list[m]] += buf[m];
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixQEq::memory_usage()
{
double bytes;
bytes = atom->nmax*nprev*2 * sizeof(double); // s_hist & t_hist
bytes += atom->nmax*11 * sizeof(double); // storage
bytes += n_cap*2 * sizeof(int); // matrix...
bytes += m_cap * sizeof(int);
bytes += m_cap * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate fictitious charge arrays
------------------------------------------------------------------------- */
void FixQEq::grow_arrays(int nmax)
{
memory->grow(s_hist,nmax,nprev,"qeq:s_hist");
memory->grow(t_hist,nmax,nprev,"qeq:t_hist");
}
/* ----------------------------------------------------------------------
copy values within fictitious charge arrays
------------------------------------------------------------------------- */
void FixQEq::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nprev; m++) {
s_hist[j][m] = s_hist[i][m];
t_hist[j][m] = t_hist[i][m];
}
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixQEq::pack_exchange(int i, double *buf)
{
for (int m = 0; m < nprev; m++) buf[m] = s_hist[i][m];
for (int m = 0; m < nprev; m++) buf[nprev+m] = t_hist[i][m];
return nprev*2;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixQEq::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < nprev; m++) s_hist[nlocal][m] = buf[m];
for (int m = 0; m < nprev; m++) t_hist[nlocal][m] = buf[nprev+m];
return nprev*2;
}
/* ---------------------------------------------------------------------- */
double FixQEq::parallel_norm( double *v, int n )
{
int i;
double my_sum, norm_sqr;
int ii;
int *ilist;
ilist = list->ilist;
my_sum = 0.0;
norm_sqr = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_sum += v[i]*v[i];
}
MPI_Allreduce( &my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
return sqrt( norm_sqr );
}
/* ---------------------------------------------------------------------- */
double FixQEq::parallel_dot( double *v1, double *v2, int n)
{
int i;
double my_dot, res;
int ii;
int *ilist;
ilist = list->ilist;
my_dot = 0.0;
res = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_dot += v1[i] * v2[i];
}
MPI_Allreduce( &my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, world );
return res;
}
/* ---------------------------------------------------------------------- */
double FixQEq::parallel_vector_acc( double *v, int n )
{
int i;
double my_acc, res;
int ii;
int *ilist;
ilist = list->ilist;
my_acc = 0.0;
res = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_acc += v[i];
}
MPI_Allreduce( &my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, world );
return res;
}
/* ---------------------------------------------------------------------- */
void FixQEq::vector_sum( double* dest, double c, double* v,
double d, double* y, int k )
{
int kk;
int *ilist;
ilist = list->ilist;
for( --k; k>=0; --k ) {
kk = ilist[k];
if (atom->mask[kk] & groupbit)
dest[kk] = c * v[kk] + d * y[kk];
}
}
/* ---------------------------------------------------------------------- */
void FixQEq::vector_add( double* dest, double c, double* v, int k )
{
int kk;
int *ilist;
ilist = list->ilist;
for( --k; k>=0; --k ) {
kk = ilist[k];
if (atom->mask[kk] & groupbit)
dest[kk] += c * v[kk];
}
}
/* ---------------------------------------------------------------------- */

129
src/QEQ/fix_qeq.h Normal file
View File

@ -0,0 +1,129 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef LMP_FIX_QEQ_H
#define LMP_FIX_QEQ_H
#include "fix.h"
#define EV_TO_KCAL_PER_MOL 14.4
#define DANGER_ZONE 0.90
#define MIN_CAP 50
#define SAFE_ZONE 1.2
#define MIN_NBRS 100
namespace LAMMPS_NS {
class FixQEq : public Fix {
public:
FixQEq(class LAMMPS *, int, char **);
~FixQEq();
int setmask();
void init_list(int,class NeighList *);
void setup_pre_force(int);
void setup_pre_force_respa(int, int);
void pre_force_respa(int, int, int);
void min_setup_pre_force(int);
void min_pre_force(int);
// derived child classes must provide these functions
virtual void init() = 0;
virtual void pre_force(int) = 0;
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
virtual int pack_reverse_comm(int, int, double *);
virtual void unpack_reverse_comm(int, int *, double *);
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
double memory_usage();
protected:
int nevery;
int n, N, m_fill;
int n_cap, nmax, m_cap;
int pack_flag;
int nlevels_respa;
class NeighList *list;
int matvecs;
double qeq_time;
double swa, swb; // lower/upper Taper cutoff radius
double Tap[8]; // Taper function
double tolerance; // tolerance for the norm of the rel residual in CG
int maxiter; // maximum number of QEq iterations
double cutoff, cutoff_sq; // neighbor cutoff
double *chi,*eta,*gamma,*zeta,*zcore; // qeq parameters
double *chizj;
double **shld;
bigint ngroup;
// fictitious charges
double *s, *t;
double **s_hist, **t_hist;
int nprev;
typedef struct{
int n, m;
int *firstnbr;
int *numnbrs;
int *jlist;
double *val;
} sparse_matrix;
sparse_matrix H;
double *Hdia_inv;
double *b_s, *b_t;
double *b_prc, *b_prm;
double *p, *q, *r, *d;
// streitz-mintmire
double alpha;
// damped dynamics
double *qf, *q1, *q2, qdamp, qstep;
void calculate_Q();
double parallel_norm(double*, int);
double parallel_dot(double*, double*, int);
double parallel_vector_acc(double*, int);
void vector_sum(double *, double, double *, double, double *,int);
void vector_add(double *, double, double *, int);
void init_storage();
void read_file(char*);
void allocate_storage();
void deallocate_storage();
void reallocate_storage();
void allocate_matrix();
void deallocate_matrix();
void reallocate_matrix();
virtual int CG(double *, double *);
virtual void sparse_matvec(sparse_matrix *, double *,double *);
};
}
#endif

292
src/QEQ/fix_qeq_dynamic.cpp Normal file
View File

@ -0,0 +1,292 @@
/* ----------------------------------------------------------------------
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: Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_qeq_dynamic.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "pair.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg)
{
qdamp = 0.10;
qstep = 0.02;
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"qdamp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/dynamic command");
qdamp = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"qstep") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/dynamic command");
qstep = atof(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/dynamic command");
}
}
/* ---------------------------------------------------------------------- */
void FixQEqDynamic::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/dynamic requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/dynamic group has no atoms");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 1;
neighbor->requests[irequest]->full = 0;
if (tolerance < 1e-4)
if (comm->me == 0)
error->warning(FLERR,"Fix qeq/dynamic tolerance may be too small"
" for damped dynamics");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixQEqDynamic::pre_force(int vflag)
{
int i,ii,iloop,inum,*ilist;
double qmass,dtq2;
double enegchkall,enegmaxall;
double *q = atom->q;
int *mask = atom->mask;
double enegchk = 0.0;
double enegtot = 0.0;
double enegmax = 0.0;
if (update->ntimestep % nevery) return;
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
if( atom->nmax > nmax ) reallocate_storage();
inum = list->inum;
ilist = list->ilist;
qmass = 0.016;
dtq2 = 0.5*qstep*qstep/qmass;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
q1[i] = q2[i] = qf[i] = 0.0;
}
for (iloop = 0; iloop < maxiter; iloop ++ ) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
q1[i] += qf[i]*dtq2 - qdamp*q1[i];
q[i] += q1[i];
}
}
pack_flag = 1;
comm->forward_comm_fix(this);
enegtot = compute_eneg();
enegtot /= ngroup;
enegchk = enegmax = 0.0;
for (ii = 0; ii < inum ; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
q2[i] = enegtot-qf[i];
enegmax = MAX(enegmax,fabs(q2[i]));
enegchk += fabs(q2[i]);
qf[i] = q2[i];
}
}
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
enegchk = enegchkall/ngroup;
MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
enegmax = enegmaxall;
if (enegchk <= tolerance && enegmax <= 100.0*tolerance) break;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
q1[i] += qf[i]*dtq2 - qdamp*q1[i];
}
}
if (comm->me == 0) {
if (iloop == maxiter) {
char str[128];
sprintf(str,"Charges did not converge at step "BIGINT_FORMAT
": %lg",enegchk,update->ntimestep);
error->warning(FLERR,str);
}
}
}
/* ---------------------------------------------------------------------- */
double FixQEqDynamic::compute_eneg()
{
int i, j, ii, jj, inum, jnum, itype, jtype, flag;
int *ilist, *jlist, *numneigh, **firstneigh;
double eneg, enegtot;
double r, rsq, delr[3], rinv;
int nlocal = atom->nlocal;
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double *q = atom->q;
double **x = atom->x;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
qf[i] = 0.0;
}
// communicating charge force to all nodes, first forward then reverse
pack_flag = 2;
comm->forward_comm_fix(this);
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
if (mask[i] & groupbit) {
qf[i] += chi[itype] + eta[itype] * q[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
delr[0] = x[i][0] - x[j][0];
delr[1] = x[i][1] - x[j][1];
delr[2] = x[i][2] - x[j][2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
if (rsq > cutoff_sq) continue;
r = sqrt(rsq);
rinv = 1.0/r;
qf[i] += q[j] * rinv;
qf[j] += q[i] * rinv;
}
}
}
comm->reverse_comm_fix(this);
// sum charge force on each node and return it
eneg = enegtot = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
eneg += qf[i];
}
MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
return enegtot;
}
/* ---------------------------------------------------------------------- */
int FixQEqDynamic::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int m;
if( pack_flag == 1 )
for(m = 0; m < n; m++) buf[m] = atom->q[list[m]];
else if( pack_flag == 2 )
for(m = 0; m < n; m++) buf[m] = qf[list[m]];
return n;
}
/* ---------------------------------------------------------------------- */
void FixQEqDynamic::unpack_forward_comm(int n, int first, double *buf)
{
int i, m;
if( pack_flag == 1)
for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
else if( pack_flag == 2)
for(m = 0, i = first; m < n; m++, i++) qf[i] = buf[m];
}
/* ---------------------------------------------------------------------- */
int FixQEqDynamic::pack_reverse_comm(int n, int first, double *buf)
{
int i, m;
for(m = 0, i = first; m < n; m++, i++) buf[m] = qf[i];
return n;
}
/* ---------------------------------------------------------------------- */
void FixQEqDynamic::unpack_reverse_comm(int n, int *list, double *buf)
{
int m;
for(m = 0; m < n; m++) qf[list[m]] += buf[m];
}
/* ---------------------------------------------------------------------- */

46
src/QEQ/fix_qeq_dynamic.h Normal file
View File

@ -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 FIX_CLASS
FixStyle(qeq/dynamic,FixQEqDynamic)
#else
#ifndef LMP_FIX_QEQ_DYNAMIC_H
#define LMP_FIX_QEQ_DYNAMIC_H
#include "fix_qeq.h"
namespace LAMMPS_NS {
class FixQEqDynamic : public FixQEq {
public:
FixQEqDynamic(class LAMMPS *, int, char **);
~FixQEqDynamic() {}
void init();
void pre_force(int);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
private:
double compute_eneg();
};
}
#endif
#endif

187
src/QEQ/fix_qeq_point.cpp Normal file
View File

@ -0,0 +1,187 @@
/* ----------------------------------------------------------------------
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: Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_qeq_point.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqPoint::FixQEqPoint(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
void FixQEqPoint::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/point requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/point group has no atoms");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 1;
neighbor->requests[irequest]->full = 0;
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixQEqPoint::pre_force(int vflag)
{
if (update->ntimestep % nevery) return;
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
if( atom->nmax > nmax ) reallocate_storage();
if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
reallocate_matrix();
init_matvec();
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
}
/* ---------------------------------------------------------------------- */
void FixQEqPoint::init_matvec()
{
compute_H();
int nn, ii, i;
int *ilist;
nn = list->inum;
ilist = list->ilist;
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
}
/* ---------------------------------------------------------------------- */
void FixQEqPoint::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int i, j, ii, jj, flag;
double **x, SMALL = 0.0001;
double dx, dy, dz, r_sqr, r;
int *type = atom->type;
tagint *tag = atom->tag;
x = atom->x;
int *mask = atom->mask;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// fill in the H matrix
m_fill = 0;
r_sqr = 0;
for( ii = 0; ii < inum; ii++ ) {
i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
dz = x[j][2] - x[i][2];
r_sqr = dx*dx + dy*dy + dz*dz;
flag = 0;
if (r_sqr <= cutoff_sq) {
if (j < n) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
else if (fabs(dz) < SMALL) {
if (dy > SMALL) flag = 1;
else if (fabs(dy) < SMALL && dx > SMALL)
flag = 1;
}
}
}
if( flag ) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);
H.val[m_fill] = 1.0/r;
m_fill++;
}
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
}
}
if (m_fill >= H.m) {
char str[128];
sprintf(str,"H matrix size has been exceeded: m_fill=%d H.m=%d\n",
m_fill, H.m );
error->warning(FLERR,str);
error->all(FLERR,"Fix qeq/point has insufficient QEq matrix size");
}
}
/* ---------------------------------------------------------------------- */

41
src/QEQ/fix_qeq_point.h Normal file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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 FIX_CLASS
FixStyle(qeq/point,FixQEqPoint)
#else
#ifndef LMP_FIX_QEQ_POINT_H
#define LMP_FIX_QEQ_POINT_H
#include "fix_qeq.h"
namespace LAMMPS_NS {
class FixQEqPoint : public FixQEq {
public:
FixQEqPoint(class LAMMPS *, int, char **);
~FixQEqPoint() {}
void init();
void pre_force(int);
private:
void init_matvec();
void compute_H();
};
}
#endif
#endif

249
src/QEQ/fix_qeq_shielded.cpp Normal file
View File

@ -0,0 +1,249 @@
/* ----------------------------------------------------------------------
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: Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_qeq_shielded.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
void FixQEqShielded::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/shielded requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/shielded group has no atoms");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 1;
neighbor->requests[irequest]->full = 0;
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding");
init_shielding();
int i;
for (i = 1; i <= ntypes; i++) {
if (gamma[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/shielded");
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixQEqShielded::init_shielding()
{
int i,j;
double d7, swa2, swa3, swb2, swb3;
int ntypes = atom->ntypes;
for( i = 1; i <= ntypes; ++i )
for( j = 1; j <= ntypes; ++j )
shld[i][j] = pow( gamma[i] * gamma[j], -1.5 );
if (fabs(swa) > 0.01 && comm->me == 0)
error->warning(FLERR,"Fix qeq has non-zero lower Taper radius cutoff");
if (swb < 0)
error->all(FLERR, "Fix qeq has negative upper Taper radius cutoff");
else if (swb < 5 && comm->me == 0)
error->warning(FLERR,"Fix qeq has very low Taper radius cutoff");
d7 = pow( swb - swa, 7 );
swa2 = swa*swa;
swa3 = swa2*swa;
swb2 = swb*swb;
swb3 = swb2*swb;
Tap[7] = 20.0 / d7;
Tap[6] = -70.0 * (swa + swb) / d7;
Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7;
Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3 ) / d7;
Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3 ) / d7;
Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7;
Tap[1] = 140.0 * swa3 * swb3 / d7;
Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 +
7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7;
}
/* ---------------------------------------------------------------------- */
void FixQEqShielded::pre_force(int vflag)
{
if (update->ntimestep % nevery) return;
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
if( atom->nmax > nmax ) reallocate_storage();
if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
reallocate_matrix();
init_matvec();
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
}
/* ---------------------------------------------------------------------- */
void FixQEqShielded::init_matvec()
{
compute_H();
int nn, ii, i;
int *ilist;
nn = list->inum;
ilist = list->ilist;
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
}
/* ---------------------------------------------------------------------- */
void FixQEqShielded::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int i, j, ii, jj, flag;
double **x, SMALL = 0.0001;
double dx, dy, dz, r_sqr, r;
int *type = atom->type;
tagint *tag = atom->tag;
x = atom->x;
int *mask = atom->mask;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// fill in the H matrix
m_fill = 0;
r_sqr = 0;
for( ii = 0; ii < inum; ii++ ) {
i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
dz = x[j][2] - x[i][2];
r_sqr = dx*dx + dy*dy + dz*dz;
flag = 0;
if (r_sqr <= cutoff_sq) {
if (j < n) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
else if (fabs(dz) < SMALL) {
if (dy > SMALL) flag = 1;
else if (fabs(dy) < SMALL && dx > SMALL)
flag = 1;
}
}
}
if( flag ) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);
H.val[m_fill] = calculate_H( r, shld[type[i]][type[j]] );
m_fill++;
}
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
}
}
if (m_fill >= H.m) {
char str[128];
sprintf(str,"H matrix size has been exceeded: m_fill=%d H.m=%d\n",
m_fill, H.m );
error->warning(FLERR,str);
error->all(FLERR,"Fix qeq/shielded has insufficient QEq matrix size");
}
}
/* ---------------------------------------------------------------------- */
double FixQEqShielded::calculate_H( double r, double gamma )
{
double Taper, denom;
Taper = Tap[7] * r + Tap[6];
Taper = Taper * r + Tap[5];
Taper = Taper * r + Tap[4];
Taper = Taper * r + Tap[3];
Taper = Taper * r + Tap[2];
Taper = Taper * r + Tap[1];
Taper = Taper * r + Tap[0];
denom = r * r * r + gamma;
denom = pow(denom,0.3333333333333);
return Taper * EV_TO_KCAL_PER_MOL / denom;
}

43
src/QEQ/fix_qeq_shielded.h Normal file
View File

@ -0,0 +1,43 @@
/* ----------------------------------------------------------------------
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 FIX_CLASS
FixStyle(qeq/shielded,FixQEqShielded)
#else
#ifndef LMP_FIX_QEQ_SHIELDED_H
#define LMP_FIX_QEQ_SHIELDED_H
#include "fix_qeq.h"
namespace LAMMPS_NS {
class FixQEqShielded : public FixQEq {
public:
FixQEqShielded(class LAMMPS *, int, char **);
~FixQEqShielded() {}
void init();
void pre_force(int);
private:
void init_shielding();
void init_matvec();
void compute_H();
double calculate_H(double,double);
};
}
#endif
#endif

431
src/QEQ/fix_qeq_slater.cpp Normal file
View File

@ -0,0 +1,431 @@
/* ----------------------------------------------------------------------
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: Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_qeq_slater.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "kspace.h"
#include "group.h"
#include "pair.h"
#include "respa.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg)
{
alpha = 0.20;
// optional arg
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"alpha") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/slater command");
alpha = atof(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/slater command");
}
}
/* ---------------------------------------------------------------------- */
void FixQEqSlater::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/slater requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/slater group has no atoms");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
int ntypes = atom->ntypes;
for (int i = 1; i <= ntypes; i++) {
if (zeta[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/slater");
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixQEqSlater::pre_force(int vflag)
{
if (update->ntimestep % nevery) return;
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
if( atom->nmax > nmax ) reallocate_storage();
if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
reallocate_matrix();
init_matvec();
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
}
/* ---------------------------------------------------------------------- */
void FixQEqSlater::init_matvec()
{
compute_H();
int nn, ii, i;
int *ilist;
nn = list->inum;
ilist = list->ilist;
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm_fix(this); //Dist_vector( s );
pack_flag = 3;
comm->forward_comm_fix(this); //Dist_vector( t );
}
/* ---------------------------------------------------------------------- */
void FixQEqSlater::compute_H()
{
int i, j, ii, jj, inum, jnum, itype, jtype, flag;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, rsq, delr[3];
double zei, zej, zj, zjtmp;
double SMALL = 0.0001;
int *tag = atom->tag;
int *type = atom->type;
double **x = atom->x;
double *q = atom->q;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
int nlocal = atom->nlocal;
m_fill = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
zei = zeta[itype];
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
zjtmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
zej = zeta[jtype];
zj = zcore[jtype];
delr[0] = x[i][0] - x[j][0];
delr[1] = x[i][1] - x[j][1];
delr[2] = x[i][2] - x[j][2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
if (rsq > cutoff_sq) continue;
r = sqrt(rsq);
H.jlist[m_fill] = j;
H.val[m_fill] = calculate_H(zei, zej, zj, r, zjtmp);
m_fill++;
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
chizj[i] = zjtmp;
}
if (m_fill >= H.m) {
char str[128];
sprintf(str,"H matrix size has been exceeded: m_fill=%d H.m=%d\n",
m_fill, H.m );
error->warning(FLERR,str);
error->all(FLERR,"Fix qeq/slater has insufficient QEq matrix size");
}
}
/* ---------------------------------------------------------------------- */
double FixQEqSlater::calculate_H(double zei, double zej, double zj,
double r, double &zjtmp)
{
double rinv = 1.0/r;
double rinv2 = rinv*rinv;
double exp2zir = exp(-2.0*zei*r);
double zei2 = zei*zei;
double zei4 = zei2*zei2;
double zei6 = zei2*zei4;
double exp2zjr = exp(-2.0*zej*r);
double zej2 = zej*zej;
double zej4 = zej2*zej2;
double zej6 = zej2*zej4;
double sm1 = 11.0/8.0;
double sm2 = 3.00/4.0;
double sm3 = 1.00/6.0;
double erfcr = erfc(alpha*r);
double qqrd2e = force->qqrd2e;
double etmp, etmp1, etmp2, etmp3, etmp4;
double e1, e2, e3, e4;
double ci_jfi, ci_fifj;
e1 = e2 = e3 = e4 = 0.0;
etmp = etmp1 = etmp2 = etmp3 = etmp4 = 0.0;
ci_jfi = -zei*exp2zir - rinv*exp2zir;
if (zei == zej) {
ci_fifj = -exp2zir*(rinv + zei*(sm1 + sm2*zei*r + sm3*zei2*r*r));
} else {
e1 = zei*zej4/((zei+zej)*(zei+zej)*(zei-zej)*(zei-zej));
e2 = zej*zei4/((zei+zej)*(zei+zej)*(zej-zei)*(zej-zei));
e3 = (3.0*zei2*zej4-zej6) /
((zei+zej)*(zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)*(zei-zej));
e4 = (3.0*zej2*zei4-zei6) /
((zei+zej)*(zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)*(zej-zei));
ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r);
}
etmp1 = 1.00 * (ci_jfi - ci_fifj);
etmp2 = 0.50 * (ci_fifj + erfcr*rinv);
zjtmp += qqrd2e * zj * etmp1;
return qqrd2e * etmp2;
}
/* ---------------------------------------------------------------------- */
double FixQEqSlater::calculate_H_wolf(double zei, double zej, double zj,
double r, double &zjtmp)
{
double rinv = 1.0/r;
double rinv2 = rinv*rinv;
double exp2zir = exp(-2.0*zei*r);
double zei2 = zei*zei;
double zei4 = zei2*zei2;
double zei6 = zei2*zei4;
double exp2zjr = exp(-2.0*zej*r);
double zej2 = zej*zej;
double zej4 = zej2*zej2;
double zej6 = zej2*zej4;
double sm1 = 11.0/8.0;
double sm2 = 3.00/4.0;
double sm3 = 1.00/6.0;
double e1, e2, e3, e4;
double rc = cutoff;
double rcinv = 1.0/rc;
double rcinv2 = rcinv*rcinv;
double exp2zirsh = exp(-2.0*zei*rc);
double exp2zjrsh = exp(-2.0*zej*rc);
double eshift, fshift, ci_jfi, ci_fifj;
double etmp1, etmp2, etmp3;
double a = alpha;
double erfcr = erfc(a*r);
double erfcrc = erfc(a*rc);
double qqrd2e = force->qqrd2e;
etmp1 = etmp2 = etmp3 = 0.0;
e1 = e2 = e3 = e4 = 0.0;
eshift = -zei*exp2zirsh - rcinv*exp2zirsh;
fshift = 2.0*zei2*exp2zirsh + rcinv2*exp2zirsh + 2.0*zei*rcinv*exp2zirsh;
ci_jfi = -zei*exp2zir - rinv*exp2zir - eshift - (r-rc)*fshift;
if (zei == zej) {
eshift = -exp2zirsh*(rcinv + zei*(sm1 + sm2*zei*rc + sm3*zei2*rc*rc));
ci_fifj = -exp2zir*(rinv + zei*(sm1 + sm2*zei*r + sm3*zei2*r*r))
- eshift - (r-rc)*fshift;
} else {
e1 = zei*zej4/((zei+zej)*(zei+zej)*(zei-zej)*(zei-zej));
e2 = zej*zei4/((zei+zej)*(zei+zej)*(zej-zei)*(zej-zei));
e3 = (3.0*zei2*zej4-zej6) /
((zei+zej)*(zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)*(zei-zej));
e4 = (3.0*zej2*zei4-zei6) /
((zei+zej)*(zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)*(zej-zei));
eshift = -exp2zirsh*(e1+e3/rc) - exp2zjrsh*(e2+e4/rc);
ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r)
- eshift - (r-rc)*fshift;
}
etmp1 = erfcr/r - erfcrc/rc;
etmp2 = 1.00 * (ci_jfi - ci_fifj);
etmp3 = 0.50 * (etmp1 + ci_fifj);
zjtmp += qqrd2e * zj * etmp2;
return qqrd2e * etmp3;
}
/* ---------------------------------------------------------------------- */
int FixQEqSlater::CG( double *b, double *x )
{
int i, j;
double tmp, alfa, beta, b_norm;
double sig_old, sig_new;
int nn, jj;
int *ilist;
nn = list->inum;
ilist = list->ilist;
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this ); //Coll_Vector( q );
vector_sum( r , 1., b, -1., q, nn );
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
d[j] = r[j] * Hdia_inv[j]; //pre-condition
}
b_norm = parallel_norm( b, nn );
sig_new = parallel_dot( r, d, nn);
for( i = 1; i < maxiter && sqrt(sig_new) / b_norm > tolerance; ++i ) {
comm->forward_comm_fix(this); //Dist_vector( d );
sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this); //Coll_vector( q );
tmp = parallel_dot( d, q, nn);
alfa = sig_new / tmp;
vector_add( x, alfa, d, nn );
vector_add( r, -alfa, q, nn );
// pre-conditioning
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
p[j] = r[j] * Hdia_inv[j];
}
sig_old = sig_new;
sig_new = parallel_dot( r, p, nn);
beta = sig_new / sig_old;
vector_sum( d, 1., p, beta, d, nn );
}
if (i >= maxiter && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq/slater CG convergence failed (%g) after %d iterations "
"at " BIGINT_FORMAT " step",sqrt(sig_new) / b_norm,i,update->ntimestep);
error->warning(FLERR,str);
}
return i;
}
/* ---------------------------------------------------------------------- */
void FixQEqSlater::sparse_matvec( sparse_matrix *A, double *x, double *b )
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
nn = atom->nlocal;
NN = atom->nlocal + atom->nghost;
ilist = list->ilist;
double r = cutoff;
double woself = 0.50*erfc(alpha*r)/r + alpha/MY_PIS;
for( i = 0; i < nn; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = (eta[atom->type[i]] - 2.0*force->qqr2e*woself) * x[i];
}
for( i = nn; i < NN; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = 0;
}
for( i = 0; i < nn; ++i ) {
if (atom->mask[i] & groupbit) {
for( itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];
}
}
}
}

44
src/QEQ/fix_qeq_slater.h Normal file
View File

@ -0,0 +1,44 @@
/* ----------------------------------------------------------------------
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 FIX_CLASS
FixStyle(qeq/slater,FixQEqSlater)
#else
#ifndef LMP_FIX_QEQ_SLATER_H
#define LMP_FIX_QEQ_SLATER_H
#include "fix_qeq.h"
namespace LAMMPS_NS {
class FixQEqSlater : public FixQEq {
public:
FixQEqSlater(class LAMMPS *, int, char **);
~FixQEqSlater() {}
void init();
void pre_force(int);
private:
void init_matvec();
int CG(double*,double*);
void sparse_matvec(sparse_matrix*,double*,double*);
void compute_H();
double calculate_H(double, double, double, double, double &);
double calculate_H_wolf(double, double, double, double, double &);
};
}
#endif
#endif