diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index 0c322634d..f822aee06 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -1,758 +1,765 @@
/* ----------------------------------------------------------------------
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 "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;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
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;
nlocal = n_cap = 0;
nall = 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;
// 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;
streitz_flag = 0;
+ qv = 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] = atom->q[i];
if (strcmp(arg[7],"coul/streitz") == 0) {
streitz_flag = 1;
} else {
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);
if (!streitz_flag) {
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::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(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");
+
+ memory->create(qv,nmax,"qeq:qv");
}
/* ---------------------------------------------------------------------- */
void FixQEq::deallocate_storage()
{
memory->destroy(s);
memory->destroy(t);
memory->destroy( Hdia_inv );
memory->destroy( b_s );
memory->destroy( b_t );
memory->destroy( p );
memory->destroy( q );
memory->destroy( r );
memory->destroy( d );
memory->destroy( chizj );
memory->destroy( qf );
memory->destroy( q1 );
memory->destroy( q2 );
+
+ memory->destroy( qv );
}
/* ---------------------------------------------------------------------- */
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;
nlocal = atom->nlocal;
n_cap = MAX( (int)(nlocal * safezone), mincap );
nall = 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)
{
if (force->newton_pair == 0)
error->all(FLERR,"QEQ with 'newton pair off' not supported");
neighbor->build_one(list);
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()
{
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
for( int i = 0; i < nall; i++ ) {
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
s[i] = t[i] = atom->q[i];
chizj[i] = 0.0;
qf[i] = 0.0;
q1[i] = 0.0;
q2[i] = 0.0;
+
+ qv[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 loop, i, ii, inum, *ilist;
double tmp, alfa, beta, b_norm;
double sig_old, sig_new;
inum = list->inum;
ilist = list->ilist;
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this );
vector_sum( r , 1., b, -1., q, inum );
for( ii = 0; ii < inum; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
d[i] = r[i] * Hdia_inv[i];
}
b_norm = parallel_norm( b, inum );
sig_new = parallel_dot( r, d, inum);
for( loop = 1; loop < maxiter && sqrt(sig_new)/b_norm > tolerance; ++loop ) {
comm->forward_comm_fix(this);
sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this);
tmp = parallel_dot( d, q, inum);
alfa = sig_new / tmp;
vector_add( x, alfa, d, inum );
vector_add( r, -alfa, q, inum );
for( ii = 0; ii < inum; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
p[i] = r[i] * Hdia_inv[i];
}
sig_old = sig_new;
sig_new = parallel_dot( r, p, inum);
beta = sig_new / sig_old;
vector_sum( d, 1., p, beta, d, inum );
}
if (loop >= 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,loop,update->ntimestep);
error->warning(FLERR,str);
}
return loop;
}
/* ---------------------------------------------------------------------- */
void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
{
int i, j, itr_j;
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
for( i = 0; i < nlocal; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = eta[ atom->type[i] ] * x[i];
}
for( i = nlocal; i < nall; ++i ) {
if (atom->mask[i] & groupbit)
b[i] = 0;
}
for( i = 0; i < nlocal; ++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, inum, ii;
int *ilist;
double u, s_sum, t_sum;
double *q = atom->q;
inum = list->inum;
ilist = list->ilist;
s_sum = parallel_vector_acc( s, inum );
t_sum = parallel_vector_acc( t, inum);
u = s_sum / t_sum;
for( ii = 0; ii < inum; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
q[i] = s[i] - u * t[i];
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 m;
}
/* ---------------------------------------------------------------------- */
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 m;
}
/* ---------------------------------------------------------------------- */
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 n, double *buf)
{
for (int m = 0; m < nprev; m++) s_hist[n][m] = buf[m];
for (int m = 0; m < nprev; m++) t_hist[n][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;
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];
}
}
/* ---------------------------------------------------------------------- */
void FixQEq::read_file(char *file)
{
int itype,ntypes;
int params_per_line = 6;
char **words = new char*[params_per_line+1];
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");
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix qeq parameter file %s",file);
error->one(FLERR,str);
}
}
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
int n,nwords,ielement,eof;
char line[MAXLINE],*ptr;
eof = ielement = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
ielement ++;
if (ielement > ntypes)
error->all(FLERR,"Invalid fix qeq parameter file");
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
itype = atoi(words[0]);
chi[itype] = atof(words[1]);
eta[itype] = atof(words[2]);
gamma[itype] = atof(words[3]);
zeta[itype] = atof(words[4]);
zcore[itype] = atof(words[5]);
}
delete [] words;
}
diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h
index 337008f9f..152ffd4e3 100644
--- a/src/QEQ/fix_qeq.h
+++ b/src/QEQ/fix_qeq.h
@@ -1,156 +1,159 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#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 nlocal, nall, 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;
int streitz_flag;
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 *p, *q, *r, *d;
// streitz-mintmire
double alpha;
// damped dynamics
double *qf, *q1, *q2, qdamp, qstep;
+ // fire
+ double *qv;
+
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
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: QEQ with 'newton pair off' not supported
See the newton command. This is a restriction to use the QEQ fixes.
W: Fix qeq CG convergence failed (%g) after %d iterations at %ld step
Self-explanatory.
E: Cannot open fix qeq parameter file %s
The specified file cannot be opened. Check that the path and name are
correct.
E: Invalid fix qeq parameter file
Element index > number of atom types.
*/
diff --git a/src/QEQ/fix_qeq_fire.cpp b/src/QEQ/fix_qeq_fire.cpp
new file mode 100644
index 000000000..7c3a9409c
--- /dev/null
+++ b/src/QEQ/fix_qeq_fire.cpp
@@ -0,0 +1,357 @@
+/* ----------------------------------------------------------------------
+ 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_fire.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 "pair_comb.h"
+#include "pair_comb3.h"
+#include "kspace.h"
+#include "respa.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define DELAYSTEP 0
+#define DT_GROW 1.1
+#define DT_SHRINK 0.5
+#define ALPHA0 0.8
+#define ALPHA_SHRINK 0.10
+#define TMAX 10.0
+
+/* ---------------------------------------------------------------------- */
+
+FixQEqFire::FixQEqFire(LAMMPS *lmp, int narg, char **arg) :
+ FixQEq(lmp, narg, arg)
+{
+ qdamp = 0.20;
+ qstep = 0.20;
+
+ int iarg = 8;
+ while (iarg < narg) {
+
+ if (strcmp(arg[iarg],"qdamp") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/fire 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/fire command");
+ qstep = atof(arg[iarg+1]);
+ iarg += 2;
+ } else error->all(FLERR,"Illegal fix qeq/fire command");
+ }
+
+ comb = NULL;
+ comb3 = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEqFire::init()
+{
+ if (!atom->q_flag)
+ error->all(FLERR,"Fix qeq/fire requires atom attribute q");
+
+ ngroup = group->count(igroup);
+ if (ngroup == 0) error->all(FLERR,"Fix qeq/fire group has no atoms");
+
+ int irequest = neighbor->request(this,instance_me);
+ 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/fire tolerance may be too small"
+ " for damped fires");
+
+ if (strstr(update->integrate_style,"respa"))
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+ comb = (PairComb *) force->pair_match("comb",1);
+ comb3 = (PairComb3 *) force->pair_match("comb3",1);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEqFire::pre_force(int vflag)
+{
+ int inum, *ilist;
+ int i,ii,iloop,loopmax;
+ int *mask = atom->mask;
+
+ double *q = atom->q;
+ double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
+ double scale1,scale2;
+ double dtvone,dtv;
+ double enegtot,enegchk,enegmax;
+ double alpha = qdamp;
+ double dt, dtmax;
+ double enegchkall,enegmaxall;
+ bigint ntimestep = update->ntimestep;
+ bigint last_negative = 0;
+
+ if (ntimestep % nevery) return;
+
+ if( atom->nmax > nmax ) reallocate_storage();
+
+ inum = list->inum;
+ ilist = list->ilist;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qv[i] = 0.0;
+ }
+
+ dt = qstep;
+ dtmax = TMAX * dt;
+
+ for (iloop = 0; iloop < maxiter; iloop ++ ) {
+ pack_flag = 1;
+ comm->forward_comm_fix(this);
+
+ if (comb) {
+ comb->yasu_char(qf,igroup);
+ enegtot = comb->enegtot / ngroup;
+ } else if (comb3) {
+ comb3->combqeq(qf,igroup);
+ enegtot = comb3->enegtot / ngroup;
+ } else {
+ enegtot = compute_eneg();
+ enegtot /= ngroup;
+ }
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qf[i] -= enegtot; // Enforce adiabatic
+ }
+
+ // FIRE minimization algorithm
+ // vdotfall = v dot f = qv dot qf
+ vdotf = 0.0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ vdotf += (qv[i]*qf[i]);
+ }
+ MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
+
+ if (vdotfall > 0.0) {
+ vdotv = fdotf = 0.0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ vdotv += qv[i]*qv[i];
+ fdotf += qf[i]*qf[i];
+ }
+ MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
+ MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
+
+ scale1 = 1.0 - alpha;
+ if (fdotfall == 0.0) scale2 = 0.0;
+ else scale2 = alpha * sqrt(vdotvall/fdotfall);
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qv[i] = scale1*qv[i] + scale2*qf[i];
+ }
+ if (ntimestep - last_negative > DELAYSTEP) {
+ dt = MIN(dt*DT_GROW,dtmax);
+ alpha *= ALPHA_SHRINK;
+ }
+ } else {
+ last_negative = ntimestep;
+ dt *= DT_SHRINK;
+ alpha = ALPHA0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qv[i] = 0.0;
+ }
+ }
+
+ // limit timestep so no charges change more than dmax
+ dtvone = dt;
+ double dmax = 0.1;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ vmax = MAX(fabs(qv[i]),0);
+ if (dtvone*vmax > dmax) dtvone = dmax/vmax;
+ }
+ MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
+ //dtv = dt;
+
+ // Euler integration step
+ enegchk = 0.0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ q[i] -= dtv * qv[i];
+ qv[i] += dtv * qf[i];
+ enegchk += fabs(qf[i]);
+ }
+ MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
+ enegchk = enegchkall / ngroup;
+
+ if (enegchk < tolerance) break;
+ }
+
+ if (comm->me == 0) {
+ if (iloop == maxiter) {
+ char str[128];
+ sprintf(str,"Charges did not converge at step "BIGINT_FORMAT
+ ": %lg",update->ntimestep,enegchk);
+ error->warning(FLERR,str);
+ }
+ }
+
+ if (force->kspace) force->kspace->qsum_qsq();
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixQEqFire::compute_eneg()
+{
+ int i, j, ii, jj, inum, jnum, itype;
+ int *ilist, *jlist, *numneigh, **firstneigh;
+ double eneg, enegtot;
+ double r, rsq, delr[3], rinv;
+
+ 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];
+ j &= NEIGHMASK;
+
+ 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;
+ }
+ }
+ }
+
+ pack_flag = 2;
+ 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 FixQEqFire::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 m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEqFire::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 FixQEqFire::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 m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEqFire::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int m;
+
+ for(m = 0; m < n; m++) qf[list[m]] += buf[m];
+}
+
+/* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_fire.h b/src/QEQ/fix_qeq_fire.h
new file mode 100644
index 000000000..9fd81176e
--- /dev/null
+++ b/src/QEQ/fix_qeq_fire.h
@@ -0,0 +1,75 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(qeq/fire,FixQEqFire)
+
+#else
+
+#ifndef LMP_FIX_QEQ_FIRE_H
+#define LMP_FIX_QEQ_FIRE_H
+
+#include "fix_qeq.h"
+
+namespace LAMMPS_NS {
+
+class FixQEqFire : public FixQEq {
+ public:
+ FixQEqFire(class LAMMPS *, int, char **);
+ ~FixQEqFire() {}
+ 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();
+
+ class PairComb *comb;
+ class PairComb3 *comb3;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix qeq/fire requires atom attribute q
+
+Self-explanatory.
+
+E: Fix qeq/fire group has no atoms
+
+Self-explanatory.
+
+W: Fix qeq/fire tolerance may be too small for damped fires
+
+Self-explanatory.
+
+W: Charges did not converge at step %ld: %lg
+
+Self-explanatory.
+
+*/