diff --git a/src/Makefile b/src/Makefile index 38a1c15fa6..d919520ffb 100755 --- a/src/Makefile +++ b/src/Makefile @@ -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 \ diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp new file mode 100644 index 0000000000..b66d31b03d --- /dev/null +++ b/src/QEQ/fix_qeq.cpp @@ -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_jfirstnbr[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]; + } +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h new file mode 100644 index 0000000000..68176fd668 --- /dev/null +++ b/src/QEQ/fix_qeq.h @@ -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 diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp new file mode 100644 index 0000000000..4f7cb0d994 --- /dev/null +++ b/src/QEQ/fix_qeq_dynamic.cpp @@ -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]; +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_dynamic.h b/src/QEQ/fix_qeq_dynamic.h new file mode 100644 index 0000000000..12481722db --- /dev/null +++ b/src/QEQ/fix_qeq_dynamic.h @@ -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 diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp new file mode 100644 index 0000000000..9dae0655c5 --- /dev/null +++ b/src/QEQ/fix_qeq_point.cpp @@ -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"); + } +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_point.h b/src/QEQ/fix_qeq_point.h new file mode 100644 index 0000000000..c73baf9040 --- /dev/null +++ b/src/QEQ/fix_qeq_point.h @@ -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 diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp new file mode 100644 index 0000000000..dfc110f75e --- /dev/null +++ b/src/QEQ/fix_qeq_shielded.cpp @@ -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; + +} diff --git a/src/QEQ/fix_qeq_shielded.h b/src/QEQ/fix_qeq_shielded.h new file mode 100644 index 0000000000..3c882bde00 --- /dev/null +++ b/src/QEQ/fix_qeq_shielded.h @@ -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 diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp new file mode 100644 index 0000000000..d13669f2d1 --- /dev/null +++ b/src/QEQ/fix_qeq_slater.cpp @@ -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_jfirstnbr[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]; + } + } + } + +} diff --git a/src/QEQ/fix_qeq_slater.h b/src/QEQ/fix_qeq_slater.h new file mode 100644 index 0000000000..167334093e --- /dev/null +++ b/src/QEQ/fix_qeq_slater.h @@ -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