diff --git a/src/Makefile b/src/Makefile index 67393412df..29e0afcd45 100755 --- a/src/Makefile +++ b/src/Makefile @@ -17,8 +17,8 @@ PACKAGE = asphere class2 colloid dipole dsmc gpu granular \ kspace manybody meam molecule opt peri poems reax replica \ shock srd xtc -PACKUSER = user-misc user-ackland user-atc user-awpmd user-cd-eam user-cg-cmm \ - user-cuda user-eff user-ewaldn user-imd user-reaxc user-smd +PACKUSER = user-misc user-atc user-awpmd user-cg-cmm \ + user-cuda user-eff user-ewaldn user-reaxc PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-ACKLAND/Install.sh b/src/USER-ACKLAND/Install.sh deleted file mode 100644 index d3cef84739..0000000000 --- a/src/USER-ACKLAND/Install.sh +++ /dev/null @@ -1,15 +0,0 @@ -# Install/unInstall package files in LAMMPS - -if (test $1 = 1) then - - cp compute_ackland_atom.cpp .. - - cp compute_ackland_atom.h .. - -elif (test $1 = 0) then - - rm -f ../compute_ackland_atom.cpp - - rm -f ../compute_ackland_atom.h - -fi diff --git a/src/USER-ACKLAND/README b/src/USER-ACKLAND/README deleted file mode 100644 index d44890ed3b..0000000000 --- a/src/USER-ACKLAND/README +++ /dev/null @@ -1,26 +0,0 @@ -The files in this directory are a user-contributed package for LAMMPS. - -The person who created these files is Gerolf Ziegenhain -(gerolf@ziegenhain.com). Contact him directly if you have questions. - -This package implements a "compute ackland/atom" command which can be -used in a LAMMPS input script. Like other per-atom compute commands, -the results can be accessed when dumping atom information to a file, -or by other fixes that do averaging of various kinds. See the -documentation files for these commands for details. - -The Ackland computation is a means of detecting local lattice -structure around an atom, as described in G. Ackland, -PRB(2006)73:054104. - -The output is a number with the following mapping: - -enum{UNKNOWN,BCC,FCC,HCP,ICO}; - -or in other words: - -0 == UNKNOWN -1 == BCC -2 == FCC -3 == HCP -4 == ICO diff --git a/src/USER-CD-EAM/Install.sh b/src/USER-CD-EAM/Install.sh deleted file mode 100755 index 614d620f58..0000000000 --- a/src/USER-CD-EAM/Install.sh +++ /dev/null @@ -1,15 +0,0 @@ -# Install/Uninstall package files in LAMMPS - -if (test $1 = 1) then - - cp pair_cdeam.cpp .. - - cp pair_cdeam.h .. - -elif (test $1 = 0) then - - rm -f ../pair_cdeam.cpp - - rm -f ../pair_cdeam.h - -fi diff --git a/src/USER-CD-EAM/README b/src/USER-CD-EAM/README deleted file mode 100644 index aadf395d5e..0000000000 --- a/src/USER-CD-EAM/README +++ /dev/null @@ -1,29 +0,0 @@ -The files in this directory are a user-contributed package for LAMMPS. - -The person who created these files is Alexander Stukowski of the -Technical University of Darmstadt, Germany Department of Materials -Science (stukowski at mm.tu-darmstadt.de). Contact him directly if -you have questions. - -This package implements the concentration-dependent EAM (CD-EAM) -potential for multi-component systems. - -Specifically, the pair style implements an improved version of the -original CD-EAM formalism. The new version (a.k.a. one-site model; -cdeamVersion==1) has been published in - - A. Stukowski, B. Sadigh, P. Erhart and A. Caro - Efficient implementation of the concentration-dependent embedded - atom method for molecular-dynamics and Monte-Carlo simulations - Model. Simul. Mater. Sci. Eng., 2009, 075005 - -This new formulation is more efficient for MD and Monte-Carlo -simulations and is the default. - -The original formulation (a.k.a. two-site model; cdeamVersion==2) is -also implemented and has been published in - - A. Caro, D. A. Crowson and M. Caro - Classical Many-Body Potential for Concentrated Alloys and the - Inversion of Order in Iron-Chromium Alloys - Phys. Rev. Lett., APS, 2005, 95, 075702 diff --git a/src/USER-CD-EAM/pair_cdeam.cpp b/src/USER-CD-EAM/pair_cdeam.cpp deleted file mode 100644 index 2d0e224f1c..0000000000 --- a/src/USER-CD-EAM/pair_cdeam.cpp +++ /dev/null @@ -1,642 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Alexander Stukowski - Technical University of Darmstadt, - Germany Department of Materials Science -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_cdeam.h" -#include "atom.h" -#include "force.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -// This is for debugging purposes. The ASSERT() macro is used in the code to check -// if everything runs as expected. Change this to #if 0 if you don't need the checking. -#if 0 - #define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop()) - - inline void my_noop() {} - inline void my_failure(Error* error, const char* file, int line) { - char str[1024]; - sprintf(str,"Assertion failure: File %s, line %i", file, line); - error->one(str); - } -#else - #define ASSERT(cond) -#endif - -#define MAXLINE 1024 // This sets the maximum line length in EAM input files. - -PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion) -{ - single_enable = 0; - rhoB = NULL; - D_values = NULL; - hcoeff = NULL; - - // Set communication buffer sizes needed by this pair style. - if(cdeamVersion == 1) { - comm_forward = 4; - comm_reverse = 3; - } - else if(cdeamVersion == 2) { - comm_forward = 3; - comm_reverse = 2; - } - else { - error->all("Invalid CD-EAM potential version."); - } -} - -PairCDEAM::~PairCDEAM() -{ - memory->destroy(rhoB); - memory->destroy(D_values); - if(hcoeff) delete[] hcoeff; -} - -void PairCDEAM::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rhoip,rhojp,recip,phi; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; - - // Grow per-atom arrays if necessary - if(atom->nmax > nmax) { - memory->destroy(rho); - memory->destroy(fp); - memory->destroy(rhoB); - memory->destroy(D_values); - nmax = atom->nmax; - memory->create(rho,nmax,"pair:rho"); - memory->create(rhoB,nmax,"pair:rhoB"); - memory->create(fp,nmax,"pair:fp"); - memory->create(D_values,nmax,"pair:D_values"); - } - - double **x = atom->x; - double **f = atom->f; - int *type = atom->type; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // Zero out per-atom arrays. - int m = nlocal + atom->nghost; - for(i = 0; i < m; i++) { - rho[i] = 0.0; - rhoB[i] = 0.0; - D_values[i] = 0.0; - } - - // Stage I - - // Compute rho and rhoB at each local atom site. - // Additionally calculate the D_i values here if we are using the one-site formulation. - // For the two-site formulation we have to calculate the D values in an extra loop (Stage II). - for(ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for(jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if(rsq < cutforcesq) { - jtype = type[j]; - double r = sqrt(rsq); - const EAMTableIndex index = radiusToTableIndex(r); - double localrho = RhoOfR(index, jtype, itype); - rho[i] += localrho; - if(jtype == speciesB) rhoB[i] += localrho; - if(newton_pair || j < nlocal) { - localrho = RhoOfR(index, itype, jtype); - rho[j] += localrho; - if(itype == speciesB) rhoB[j] += localrho; - } - - if(cdeamVersion == 1 && itype != jtype) { - // Note: if the i-j interaction is not concentration dependent (because either - // i or j are not species A or B) then its contribution to D_i and D_j should - // be ignored. - // This if-clause is only required for a ternary. - if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) { - double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r); - D_values[i] += Phi_AB; - if(newton_pair || j < nlocal) - D_values[j] += Phi_AB; - } - } - } - } - } - - // Communicate and sum densities. - if(newton_pair) { - communicationStage = 1; - comm->reverse_comm_pair(this); - } - - // fp = derivative of embedding energy at each atom - // phi = embedding energy at each atom - for(ii = 0; ii < inum; ii++) { - i = ilist[ii]; - EAMTableIndex index = rhoToTableIndex(rho[i]); - fp[i] = FPrimeOfRho(index, type[i]); - if(eflag) { - phi = FofRho(index, type[i]); - if (eflag_global) eng_vdwl += phi; - if (eflag_atom) eatom[i] += phi; - } - } - - // Communicate derivative of embedding function and densities - // and D_values (this for one-site formulation only). - communicationStage = 2; - comm->forward_comm_pair(this); - - // The electron densities may not drop to zero because then the concentration would no longer be defined. - // But the concentration is not needed anyway if there is no interaction with another atom, which is the case - // if the electron density is exactly zero. That's why the following lines have been commented out. - // - //for(i = 0; i < nlocal + atom->nghost; i++) { - // if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB)) - // error->one("CD-EAM potential routine: Detected atom with zero electron density."); - //} - - // Stage II - // This is only required for the original two-site formulation of the CD-EAM potential. - - if(cdeamVersion == 2) { - // Compute intermediate value D_i for each atom. - for(ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - // This code line is required for ternary alloys. - if(itype != speciesA && itype != speciesB) continue; - - double x_i = rhoB[i] / rho[i]; // Concentration at atom i. - - for(jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - jtype = type[j]; - if(itype == jtype) continue; - - // This code line is required for ternary alloys. - if(jtype != speciesA && jtype != speciesB) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if(rsq < cutforcesq) { - double r = sqrt(rsq); - const EAMTableIndex index = radiusToTableIndex(r); - - // The concentration independent part of the cross pair potential. - double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r); - - // Average concentration of two sites - double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]); - - // Calculate derivative of h(x_ij) polynomial function. - double h_prime = evalHprime(x_ij); - - D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]); - if(newton_pair || j < nlocal) - D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]); - } - } - } - - // Communicate and sum D values. - if(newton_pair) { - communicationStage = 3; - comm->reverse_comm_pair(this); - } - communicationStage = 4; - comm->forward_comm_pair(this); - } - - // Stage III - - // Compute force acting on each atom. - for(ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - - jlist = firstneigh[i]; - jnum = numneigh[i]; - - // Concentration at site i - double x_i = -1.0; // The value -1 indicates: no concentration dependence for all interactions of atom i. - // It will be replaced by the concentration at site i if atom i is either A or B. - - double D_i, h_prime_i; - - // This if-clause is only required for ternary alloys. - if((itype == speciesA || itype == speciesB) && rho[i] != 0.0) { - - // Compute local concentration at site i. - x_i = rhoB[i]/rho[i]; - ASSERT(x_i >= 0 && x_i<=1.0); - - if(cdeamVersion == 1) { - // Calculate derivative of h(x_i) polynomial function. - h_prime_i = evalHprime(x_i); - D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]); - } - else if(cdeamVersion == 2) { - D_i = D_values[i]; - } - else ASSERT(false); - } - - for(jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if(rsq < cutforcesq) { - jtype = type[j]; - double r = sqrt(rsq); - const EAMTableIndex index = radiusToTableIndex(r); - - // rhoip = derivative of (density at atom j due to atom i) - // rhojp = derivative of (density at atom i due to atom j) - // psip needs both fp[i] and fp[j] terms since r_ij appears in two - // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) - // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip - rhoip = RhoPrimeOfR(index, itype, jtype); - rhojp = RhoPrimeOfR(index, jtype, itype); - fpair = fp[i]*rhojp + fp[j]*rhoip; - recip = 1.0/r; - - double x_j = -1; // The value -1 indicates: no concentration dependence for this i-j pair - // because atom j is not of species A nor B. - - // This code line is required for ternary alloy. - if(jtype == speciesA || jtype == speciesB) { - ASSERT(rho[i] != 0.0); - ASSERT(rho[j] != 0.0); - - // Compute local concentration at site j. - x_j = rhoB[j]/rho[j]; - ASSERT(x_j >= 0 && x_j<=1.0); - - double D_j; - if(cdeamVersion == 1) { - // Calculate derivative of h(x_j) polynomial function. - double h_prime_j = evalHprime(x_j); - D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]); - } - else if(cdeamVersion == 2) { - D_j = D_values[j]; - } - else ASSERT(false); - - double t2 = -rhoB[j]; - if(itype == speciesB) t2 += rho[j]; - fpair += D_j * rhoip * t2; - } - - // This if-clause is only required for a ternary alloy. - // Actually we don't need it at all because D_i should be zero anyway if - // atom i has no concentration dependent interactions (because it is not species A or B). - if(x_i != -1.0) { - double t1 = -rhoB[i]; - if(jtype == speciesB) t1 += rho[i]; - fpair += D_i * rhojp * t1; - } - - double phip; - double phi = PhiOfR(index, itype, jtype, recip, phip); - if(itype == jtype || x_i == -1.0 || x_j == -1.0) { - // Case of no concentration dependence. - fpair += phip; - } - else { - // We have a concentration dependence for the i-j interaction. - double h; - if(cdeamVersion == 1) { - // Calculate h(x_i) polynomial function. - double h_i = evalH(x_i); - // Calculate h(x_j) polynomial function. - double h_j = evalH(x_j); - h = 0.5 * (h_i + h_j); - } - else if(cdeamVersion == 2) { - // Average concentration. - double x_ij = 0.5 * (x_i + x_j); - // Calculate h(x_ij) polynomial function. - h = evalH(x_ij); - } - else ASSERT(false); - fpair += h * phip; - phi *= h; - } - - // Divide by r_ij and negate to get forces from gradient. - fpair /= -r; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if(newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if(eflag) evdwl = phi; - if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); - } - } - } - - if(vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairCDEAM::coeff(int narg, char **arg) -{ - PairEAMAlloy::coeff(narg, arg); - - // Make sure the EAM file is a CD-EAM binary alloy. - if(setfl->nelements < 2) - error->all("The EAM file must contain at least 2 elements to be used with the eam/cd pair style."); - - // Read in the coefficients of the h polynomial from the end of the EAM file. - read_h_coeff(arg[2]); - - // Determine which atom type is the A species and which is the B species in the alloy. - // By default take the first element (index 0) in the EAM file as the A species - // and the second element (index 1) in the EAM file as the B species. - speciesA = -1; - speciesB = -1; - for(int i = 1; i <= atom->ntypes; i++) { - if(map[i] == 0) { - if(speciesA >= 0) - error->all("The first element from the EAM file may only be mapped to a single atom type."); - speciesA = i; - } - if(map[i] == 1) { - if(speciesB >= 0) - error->all("The second element from the EAM file may only be mapped to a single atom type."); - speciesB = i; - } - } - if(speciesA < 0) - error->all("The first element from the EAM file must be mapped to exactly one atom type."); - if(speciesB < 0) - error->all("The second element from the EAM file must be mapped to exactly one atom type."); -} - -/* ---------------------------------------------------------------------- - Reads in the h(x) polynomial coefficients -------------------------------------------------------------------------- */ -void PairCDEAM::read_h_coeff(char *filename) -{ - if(comm->me == 0) { - // Open potential file - FILE *fp; - char line[MAXLINE]; - char nextline[MAXLINE]; - fp = fopen(filename,"r"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open EAM potential file %s", filename); - error->one(str); - } - - // h coefficients are stored at the end of the file. - // Skip to last line of file. - while(fgets(nextline, MAXLINE, fp) != NULL) { - strcpy(line, nextline); - } - char* ptr = strtok(line, " \t\n\r\f"); - int degree = atoi(ptr); - nhcoeff = degree+1; - hcoeff = new double[nhcoeff]; - int i = 0; - while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) { - hcoeff[i++] = atof(ptr); - } - if(i != nhcoeff || nhcoeff < 1) - error->one("Failed to read h(x) function coefficients from EAM file."); - - // Close the potential file. - fclose(fp); - } - - MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world); - if(comm->me != 0) hcoeff = new double[nhcoeff]; - MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world); -} - - -/* ---------------------------------------------------------------------- */ - -int PairCDEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) -{ - int i,j,m; - - m = 0; - if(communicationStage == 2) { - if(cdeamVersion == 1) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = fp[j]; - buf[m++] = rho[j]; - buf[m++] = rhoB[j]; - buf[m++] = D_values[j]; - } - return 4; - } - else if(cdeamVersion == 2) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = fp[j]; - buf[m++] = rho[j]; - buf[m++] = rhoB[j]; - } - return 3; - } - else { ASSERT(false); return 0; } - } - else if(communicationStage == 4) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = D_values[j]; - } - return 1; - } - else return 0; -} - -/* ---------------------------------------------------------------------- */ - -void PairCDEAM::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - if(communicationStage == 2) { - if(cdeamVersion == 1) { - for(i = first; i < last; i++) { - fp[i] = buf[m++]; - rho[i] = buf[m++]; - rhoB[i] = buf[m++]; - D_values[i] = buf[m++]; - } - } - else if(cdeamVersion == 2) { - for(i = first; i < last; i++) { - fp[i] = buf[m++]; - rho[i] = buf[m++]; - rhoB[i] = buf[m++]; - } - } - else ASSERT(false); - } - else if(communicationStage == 4) { - for(i = first; i < last; i++) { - D_values[i] = buf[m++]; - } - } -} - -/* ---------------------------------------------------------------------- */ -int PairCDEAM::pack_reverse_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - - if(communicationStage == 1) { - if(cdeamVersion == 1) { - for(i = first; i < last; i++) { - buf[m++] = rho[i]; - buf[m++] = rhoB[i]; - buf[m++] = D_values[i]; - } - return 3; - } - else if(cdeamVersion == 2) { - for(i = first; i < last; i++) { - buf[m++] = rho[i]; - buf[m++] = rhoB[i]; - } - return 2; - } - else { ASSERT(false); return 0; } - } - else if(communicationStage == 3) { - for(i = first; i < last; i++) { - buf[m++] = D_values[i]; - } - return 1; - } - else return 0; -} - -/* ---------------------------------------------------------------------- */ - -void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - if(communicationStage == 1) { - if(cdeamVersion == 1) { - for(i = 0; i < n; i++) { - j = list[i]; - rho[j] += buf[m++]; - rhoB[j] += buf[m++]; - D_values[j] += buf[m++]; - } - } - else if(cdeamVersion == 2) { - for(i = 0; i < n; i++) { - j = list[i]; - rho[j] += buf[m++]; - rhoB[j] += buf[m++]; - } - } - else ASSERT(false); - } - else if(communicationStage == 3) { - for(i = 0; i < n; i++) { - j = list[i]; - D_values[j] += buf[m++]; - } - } -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ -double PairCDEAM::memory_usage() -{ - double bytes = 2 * nmax * sizeof(double); - return PairEAMAlloy::memory_usage() + bytes; -} diff --git a/src/USER-CD-EAM/pair_cdeam.h b/src/USER-CD-EAM/pair_cdeam.h deleted file mode 100644 index e39923d568..0000000000 --- a/src/USER-CD-EAM/pair_cdeam.h +++ /dev/null @@ -1,230 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(eam/cd,PairCDEAM_OneSite) -PairStyle(eam/cd/old,PairCDEAM_TwoSite) - -#else - -#ifndef LMP_PAIR_CDEAM_H -#define LMP_PAIR_CDEAM_H - -#include "pair_eam_alloy.h" - -namespace LAMMPS_NS { - -class PairCDEAM : public PairEAMAlloy -{ -public: - /// Constructor. - PairCDEAM(class LAMMPS*, int cdeamVersion); - - /// Destructor. - virtual ~PairCDEAM(); - - /// Calculates the energies and forces for all atoms in the system. - void compute(int, int); - - /// Parses the pair_coeff command parameters for this pair style. - void coeff(int, char **); - - /// This is for MPI communication with neighbor nodes. - int pack_comm(int, int *, double *, int, int *); - void unpack_comm(int, int, double *); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); - - /// Reports the memory usage of this pair style to LAMMPS. - double memory_usage(); - - /// Parses the coefficients of the h polynomial from the end of the EAM file. - void read_h_coeff(char* filename); - - public: - // The public interface exposed by this potential class. - - // Evaluates the h(x) polynomial for a given local concentration x. - inline double evalH(double x) const { - double v = 0.0; - for(int i = nhcoeff-1; i >= 1; i--) { - v = (v + hcoeff[i]) * x; - } - return v + hcoeff[0]; - } - - // Calculates the derivative of the h(x) polynomial. - inline double evalHprime(double x) const { - double v = 0.0; - for(int i = nhcoeff-1; i >= 2; i--) { - v = (v + (double)i * hcoeff[i]) * x; - } - return v + hcoeff[1]; - } - - // We have two versions of the CD-EAM potential. The user can choose which one he wants to use: - // - // Version 1 - One-site concentration: The h(x_i) function depends only on the concentration at the atomic site i. - // This is a new version with a slight modification to the formula. It happens to be computationally more efficient. - // It has been published in the MSMSE 2009 paper. - // - // Version 2 - Two-site concentration: The h(x_ij) function depends on the concentrations at two atomic sites i and j. - // This is the original version from the 2005 PRL paper. - int cdeamVersion; - - // Per-atom arrays - - // The partial density of B atoms at each atom site. - double *rhoB; - - // The intermediate value D_i for each atom. - // The meaning of these values depend on the version of the CD-EAM potential used: - // - // For the one-site concentration version these are the v_i values defined in equation (21) - // of the MSMSE paper. - // - // For the old two-site concentration version these are the M_i values defined in equation (12) - // of the MSMSE paper. - double *D_values; - - // The atom type index that is considered to be the A species in the alloy. - int speciesA; - - // The atom type index that is considered to be the B species in the alloy. - int speciesB; - - protected: - - // Evaluation functions: - - // This structure specifies an entry in one of the EAM spline tables - // and the corresponding floating point part. - typedef struct { - int m; - double p; - } EAMTableIndex; - - // Converts a radius value to an index value to be used in a spline table lookup. - inline EAMTableIndex radiusToTableIndex(double r) const { - EAMTableIndex index; - index.p = r*rdr + 1.0; - index.m = static_cast(index.p); - index.m = index.m <= (nr-1) ? index.m : (nr-1); - index.p -= index.m; - index.p = index.p <= 1.0 ? index.p : 1.0; - return index; - } - - // Converts a density value to an index value to be used in a spline table lookup. - inline EAMTableIndex rhoToTableIndex(double rho) const { - EAMTableIndex index; - index.p = rho*rdrho + 1.0; - index.m = static_cast(index.p); - index.m = index.m <= (nrho-1) ? index.m : (nrho-1); - index.p -= index.m; - index.p = index.p <= 1.0 ? index.p : 1.0; - return index; - } - - // Computes the derivative of rho(r) - inline double RhoPrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { - const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; - return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; - } - - // Computes rho(r) - inline double RhoOfR(const EAMTableIndex& index, int itype, int jtype) const { - const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; - return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; - } - - // Computes the derivative of F(rho) - inline double FPrimeOfRho(const EAMTableIndex& index, int itype) const { - const double* coeff = frho_spline[type2frho[itype]][index.m]; - return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; - } - - // Computes F(rho) - inline double FofRho(const EAMTableIndex& index, int itype) const { - const double* coeff = frho_spline[type2frho[itype]][index.m]; - return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; - } - - // Computes the derivative of z2(r) - inline double Z2PrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { - const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; - return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; - } - - // Computes z2(r) - inline double Z2OfR(const EAMTableIndex& index, int itype, int jtype) const { - const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; - return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; - } - - // Computes pair potential V_ij(r). - inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR) const { - // phi = pair potential energy - // z2 = phi * r - const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; - const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; - return z2 * oneOverR; - } - - // Computes pair potential V_ij(r) and its derivative. - inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR, double& phid) const { - // phi = pair potential energy - // phip = phi' - // z2 = phi * r - // z2p = (phi * r)' = (phi' r) + phi - const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; - const double z2p = (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; - const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; - const double phi = z2 * oneOverR; - phid = z2p * oneOverR - phi * oneOverR; - return phi; - } - - // Parameters - - // h() polynomial function coefficients - double* hcoeff; - // The number of coefficients in the polynomial. - int nhcoeff; - - // This specifies the calculation stage to let the pack/unpack communication routines know - // which data to exchange. - int communicationStage; -}; - -/// The one-site concentration formulation of CD-EAM. - class PairCDEAM_OneSite : public PairCDEAM - { - public: - /// Constructor. - PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {} - }; - - /// The two-site concentration formulation of CD-EAM. - class PairCDEAM_TwoSite : public PairCDEAM - { - public: - /// Constructor. - PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {} - }; - -}; - -#endif -#endif diff --git a/src/USER-IMD/Install.sh b/src/USER-IMD/Install.sh deleted file mode 100644 index 4f009e5638..0000000000 --- a/src/USER-IMD/Install.sh +++ /dev/null @@ -1,15 +0,0 @@ -# Install/unInstall package files in LAMMPS - -if (test $1 = 1) then - - cp -p fix_imd.cpp .. - - cp -p fix_imd.h .. - -elif (test $1 = 0) then - - rm -f ../fix_imd.cpp - - rm -f ../fix_imd.h - -fi diff --git a/src/USER-IMD/README b/src/USER-IMD/README deleted file mode 100644 index 8e3238ca81..0000000000 --- a/src/USER-IMD/README +++ /dev/null @@ -1,25 +0,0 @@ -The files in this directory are a user-contributed package for LAMMPS. - -The person who created these files is Axel Kohlmeyer -(akohlmey@gmail.com). Contact him directly if you -have questions or for example scripts that use it. - -This package implements a "fix imd" command which can be used in a -LAMMPS input script. IMD stands for interactive molecular dynamics, -and allows realtime visualization and manipulation of MD simulations -through the IMD protocol, initially implemented in VMD and NAMD. - -If LAMMPS is compiled with the preprocessor flag -DLAMMPS_ASYNC_IMD -then fix imd will use posix threads to spawn a thread on MPI rank 0 -in order to offload data reading and writing from the main execution -thread and potentiall lower the inferred latencies for slow -communication links. This feature has only been tested under linux. - -There are example scripts for using this package with LAMMPS in -examples/USER/imd. Additional examples and a driver for use with the -Novint Falcon game controller as haptic device can be found at: -http://sites.google.com/site/akohlmey/software/vrpn-icms - -This software includes code developed by the Theoretical and Computational -Biophysics Group in the Beckman Institute for Advanced Science and -Technology at the University of Illinois at Urbana-Champaign. diff --git a/src/USER-IMD/fix_imd.cpp b/src/USER-IMD/fix_imd.cpp deleted file mode 100644 index 3b34ac69b6..0000000000 --- a/src/USER-IMD/fix_imd.cpp +++ /dev/null @@ -1,1476 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - The FixIMD class contains code from VMD and NAMD which is copyrighted - by the Board of Trustees of the University of Illinois and is free to - use with LAMMPS according to point 2 of the UIUC license (10% clause): - -" Licensee may, at its own expense, create and freely distribute -complimentary works that interoperate with the Software, directing others to -the TCBG server to license and obtain the Software itself. Licensee may, at -its own expense, modify the Software to make derivative works. Except as -explicitly provided below, this License shall apply to any derivative work -as it does to the original Software distributed by Illinois. Any derivative -work should be clearly marked and renamed to notify users that it is a -modified version and not the original Software distributed by Illinois. -Licensee agrees to reproduce the copyright notice and other proprietary -markings on any derivative work and to include in the documentation of such -work the acknowledgement: - - "This software includes code developed by the Theoretical and Computational - Biophysics Group in the Beckman Institute for Advanced Science and - Technology at the University of Illinois at Urbana-Champaign." - -Licensee may redistribute without restriction works with up to 1/2 of their -non-comment source code derived from at most 1/10 of the non-comment source -code developed by Illinois and contained in the Software, provided that the -above directions for notice and acknowledgement are observed. Any other -distribution of the Software or any derivative work requires a separate -license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to -negotiate an appropriate license for such distribution." -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (TempleU) - IMD API, hash, and socket code written by: John E. Stone, - Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) -------------------------------------------------------------------------- */ - -#include -#include -#include -#include - -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) -#include -#else -#include -#include -#include -#include -#include -#include -#include -#include -#include -#endif - -#include - -#include "fix_imd.h" -#include "atom.h" -#include "comm.h" -#include "update.h" -#include "respa.h" -#include "domain.h" -#include "error.h" -#include "group.h" -#include "memory.h" - -/* re-usable integer hash table code with static linkage. */ - -/** hash table top level data structure */ -typedef struct inthash_t { - struct inthash_node_t **bucket; /* array of hash nodes */ - int size; /* size of the array */ - int entries; /* number of entries in table */ - int downshift; /* shift cound, used in hash function */ - int mask; /* used to select bits for hashing */ -} inthash_t; - -/** hash table node data structure */ -typedef struct inthash_node_t { - int data; /* data in hash node */ - int key; /* key for hash lookup */ - struct inthash_node_t *next; /* next node in hash chain */ -} inthash_node_t; - -#define HASH_FAIL -1 -#define HASH_LIMIT 0.5 - -/* initialize new hash table */ -static void inthash_init(inthash_t *tptr, int buckets); -/* lookup entry in hash table */ -static int inthash_lookup(const inthash_t *tptr, int key); -/* generate list of keys for reverse lookups. */ -static int *inthash_keys(inthash_t *tptr); -/* insert an entry into hash table. */ -static int inthash_insert(inthash_t *tptr, int key, int data); -/* delete the hash table */ -static void inthash_destroy(inthash_t *tptr); -/* adapted sort for in-place sorting of map indices. */ -static void id_sort(int *idmap, int left, int right); - -/************************************************************************ - * integer hash code: - ************************************************************************/ - -/* inthash() - Hash function returns a hash number for a given key. - * tptr: Pointer to a hash table, key: The key to create a hash number for */ -static int inthash(const inthash_t *tptr, int key) { - int hashvalue; - - hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask); - if (hashvalue < 0) { - hashvalue = 0; - } - - return hashvalue; -} - -/* - * rebuild_table_int() - Create new hash table when old one fills up. - * - * tptr: Pointer to a hash table - */ -static void rebuild_table_int(inthash_t *tptr) { - inthash_node_t **old_bucket, *old_hash, *tmp; - int old_size, h, i; - - old_bucket=tptr->bucket; - old_size=tptr->size; - - /* create a new table and rehash old buckets */ - inthash_init(tptr, old_size<<1); - for (i=0; inext; - h=inthash(tptr, tmp->key); - tmp->next=tptr->bucket[h]; - tptr->bucket[h]=tmp; - tptr->entries++; - } /* while */ - } /* for */ - - /* free memory used by old table */ - free(old_bucket); - - return; -} - -/* - * inthash_init() - Initialize a new hash table. - * - * tptr: Pointer to the hash table to initialize - * buckets: The number of initial buckets to create - */ -void inthash_init(inthash_t *tptr, int buckets) { - - /* make sure we allocate something */ - if (buckets==0) - buckets=16; - - /* initialize the table */ - tptr->entries=0; - tptr->size=2; - tptr->mask=1; - tptr->downshift=29; - - /* ensure buckets is a power of 2 */ - while (tptr->sizesize<<=1; - tptr->mask=(tptr->mask<<1)+1; - tptr->downshift--; - } /* while */ - - /* allocate memory for table */ - tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *)); - - return; -} - -/* - * inthash_lookup() - Lookup an entry in the hash table and return a pointer to - * it or HASH_FAIL if it wasn't found. - * - * tptr: Pointer to the hash table - * key: The key to lookup - */ -int inthash_lookup(const inthash_t *tptr, int key) { - int h; - inthash_node_t *node; - - - /* find the entry in the hash table */ - h=inthash(tptr, key); - for (node=tptr->bucket[h]; node!=NULL; node=node->next) { - if (node->key == key) - break; - } - - /* return the entry if it exists, or HASH_FAIL */ - return(node ? node->data : HASH_FAIL); -} - - -/* - * inthash_keys() - Return a list of keys. - * NOTE: the returned list must be freed with free(3). - */ -int *inthash_keys(inthash_t *tptr) { - - int *keys; - inthash_node_t *node; - - keys = (int *)calloc(tptr->entries, sizeof(int)); - - for (int i=0; i < tptr->size; ++i) { - for (node=tptr->bucket[i]; node != NULL; node=node->next) { - keys[node->data] = node->key; - } - } - - return keys; -} - -/* - * inthash_insert() - Insert an entry into the hash table. If the entry already - * exists return a pointer to it, otherwise return HASH_FAIL. - * - * tptr: A pointer to the hash table - * key: The key to insert into the hash table - * data: A pointer to the data to insert into the hash table - */ -int inthash_insert(inthash_t *tptr, int key, int data) { - int tmp; - inthash_node_t *node; - int h; - - /* check to see if the entry exists */ - if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL) - return(tmp); - - /* expand the table if needed */ - while (tptr->entries>=HASH_LIMIT*tptr->size) - rebuild_table_int(tptr); - - /* insert the new entry */ - h=inthash(tptr, key); - node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t)); - node->data=data; - node->key=key; - node->next=tptr->bucket[h]; - tptr->bucket[h]=node; - tptr->entries++; - - return HASH_FAIL; -} - -/* - * inthash_destroy() - Delete the entire table, and all remaining entries. - * - */ -void inthash_destroy(inthash_t *tptr) { - inthash_node_t *node, *last; - int i; - - for (i=0; isize; i++) { - node = tptr->bucket[i]; - while (node != NULL) { - last = node; - node = node->next; - free(last); - } - } - - /* free the entire array of buckets */ - if (tptr->bucket != NULL) { - free(tptr->bucket); - memset(tptr, 0, sizeof(inthash_t)); - } -} - -/************************************************************************ - * integer list sort code: - ************************************************************************/ - -/* sort for integer map. initial call id_sort(idmap, 0, natoms - 1); */ -static void id_sort(int *idmap, int left, int right) -{ - int pivot, l_hold, r_hold; - - l_hold = left; - r_hold = right; - pivot = idmap[left]; - - while (left < right) { - while ((idmap[right] >= pivot) && (left < right)) - right--; - if (left != right) { - idmap[left] = idmap[right]; - left++; - } - while ((idmap[left] <= pivot) && (left < right)) - left++; - if (left != right) { - idmap[right] = idmap[left]; - right--; - } - } - idmap[left] = pivot; - pivot = left; - left = l_hold; - right = r_hold; - - if (left < pivot) - id_sort(idmap, left, pivot-1); - if (right > pivot) - id_sort(idmap, pivot+1, right); -} - -/********** API definitions of the VMD/NAMD code ************************ - * This code was taken and adapted from VMD-1.8.7/NAMD-2.7 in Sep 2009. * - * If there are any bugs or problems, please contact akohlmey@gmail.com * - ************************************************************************/ - -/*************************************************************************** - *cr - *cr (C) Copyright 1995-2009 The Board of Trustees of the - *cr University of Illinois - *cr All Rights Reserved - *cr - ***************************************************************************/ - -/* part 1: Interactive MD (IMD) API */ - -#include - -#if ( INT_MAX == 2147483647 ) -typedef int int32; -#else -typedef short int32; -#endif - -typedef struct { - int32 type; - int32 length; -} IMDheader; - -#define IMDHEADERSIZE 8 -#define IMDVERSION 2 - -typedef enum IMDType_t { - IMD_DISCONNECT, /**< close IMD connection, leaving sim running */ - IMD_ENERGIES, /**< energy data block */ - IMD_FCOORDS, /**< atom coordinates */ - IMD_GO, /**< start the simulation */ - IMD_HANDSHAKE, /**< endianism and version check message */ - IMD_KILL, /**< kill the simulation job, shutdown IMD */ - IMD_MDCOMM, /**< MDComm style force data */ - IMD_PAUSE, /**< pause the running simulation */ - IMD_TRATE, /**< set IMD update transmission rate */ - IMD_IOERROR /**< indicate an I/O error */ -} IMDType; /**< IMD command message type enumerations */ - -typedef struct { - int32 tstep; /**< integer timestep index */ - float T; /**< Temperature in degrees Kelvin */ - float Etot; /**< Total energy, in Kcal/mol */ - float Epot; /**< Potential energy, in Kcal/mol */ - float Evdw; /**< Van der Waals energy, in Kcal/mol */ - float Eelec; /**< Electrostatic energy, in Kcal/mol */ - float Ebond; /**< Bond energy, Kcal/mol */ - float Eangle; /**< Angle energy, Kcal/mol */ - float Edihe; /**< Dihedral energy, Kcal/mol */ - float Eimpr; /**< Improper energy, Kcal/mol */ -} IMDEnergies; /**< IMD simulation energy report structure */ - -/** Send control messages - these consist of a header with no subsequent data */ -static int imd_handshake(void *); /**< check endianness, version compat */ -/** Receive header and data */ -static IMDType imd_recv_header(void *, int32 *); -/** Receive MDComm-style forces, units are Kcal/mol/angstrom */ -static int imd_recv_mdcomm(void *, int32, int32 *, float *); -/** Receive energies */ -static int imd_recv_energies(void *, IMDEnergies *); -/** Receive atom coordinates. */ -static int imd_recv_fcoords(void *, int32, float *); -/** Prepare IMD data packet header */ -static void imd_fill_header(IMDheader *header, IMDType type, int32 length); -/** Write data to socket */ -static int32 imd_writen(void *s, const char *ptr, int32 n); - -/* part 2: abstracts platform-dependent routines/APIs for using sockets */ - -typedef struct { - struct sockaddr_in addr; /* address of socket provided by bind() */ - int addrlen; /* size of the addr struct */ - int sd; /* socket file descriptor */ -} imdsocket; - -static int imdsock_init(void); -static void *imdsock_create(void); -static int imdsock_bind(void *, int); -static int imdsock_listen(void *); -static void *imdsock_accept(void *); /* return new socket */ -static int imdsock_write(void *, const void *, int); -static int imdsock_read(void *, void *, int); -static int imdsock_selread(void *, int); -static int imdsock_selwrite(void *, int); -static void imdsock_shutdown(void *); -static void imdsock_destroy(void *); - -/*************************************************************** - * End of API definitions of the VMD/NAMD code. * - * The implementation follows at the end of the file. * - ***************************************************************/ - -using namespace LAMMPS_NS; - -/* struct for packed data communication of coordinates and forces. */ -struct commdata { - int tag; - float x,y,z; -}; - -/*************************************************************** - * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID imd [unwrap (on|off)] [fscale ] - ***************************************************************/ -FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 4) - error->all("Illegal fix imd command"); - - imd_port = atoi(arg[3]); - if (imd_port < 1024) - error->all("Illegal fix imd parameter: port < 1024"); - - /* default values for optional flags */ - unwrap_flag = 0; - nowait_flag = 0; - connect_msg = 1; - imd_fscale = 1.0; - imd_trate = 1; - - /* parse optional arguments */ - int argsdone = 4; - while (argsdone+1 < narg) { - if (0 == strcmp(arg[argsdone], "unwrap")) { - if (0 == strcmp(arg[argsdone+1], "on")) { - unwrap_flag = 1; - } else { - unwrap_flag = 0; - } - } else if (0 == strcmp(arg[argsdone], "nowait")) { - if (0 == strcmp(arg[argsdone+1], "on")) { - nowait_flag = 1; - } else { - nowait_flag = 0; - } - } else if (0 == strcmp(arg[argsdone], "fscale")) { - imd_fscale = atof(arg[argsdone+1]); - } else if (0 == strcmp(arg[argsdone], "trate")) { - imd_trate = atoi(arg[argsdone+1]); - } else { - error->all("Unknown fix imd parameter"); - } - ++argsdone; ++argsdone; - } - - /* sanity check on parameters */ - if (imd_trate < 1) - error->all("Illegal fix imd parameter. trate < 1."); - - bigint n = group->count(igroup); - if (n > MAXSMALLINT) error->all("Too many atoms for fix imd"); - num_coords = static_cast (n); - - MPI_Comm_rank(world,&me); - - /* initialize various imd state variables. */ - clientsock = NULL; - localsock = NULL; - nlevels_respa = 0; - imd_inactive = 0; - imd_terminate = 0; - imd_forces = 0; - force_buf = NULL; - maxbuf = 0; - msgdata = NULL; - msglen = 0; - comm_buf = NULL; - idmap = NULL; - rev_idmap = NULL; - - if (me == 0) { - /* set up incoming socket on MPI rank 0. */ - imdsock_init(); - localsock = imdsock_create(); - clientsock = NULL; - if (imdsock_bind(localsock,imd_port)) { - perror("bind to socket failed"); - imdsock_destroy(localsock); - imd_terminate = 1; - } else { - imdsock_listen(localsock); - } - } - MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); - if (imd_terminate) - error->all("LAMMPS Terminated on error in IMD."); - - /* storage required to communicate a single coordinate or force. */ - size_one = sizeof(struct commdata); - -#if defined(LAMMPS_ASYNC_IMD) - /* set up for i/o worker thread on MPI rank 0.*/ - if (me == 0) { - if (screen) - fputs("Using fix imd with asynchronous I/O.\n",screen); - if (logfile) - fputs("Using fix imd with asynchronous I/O.\n",logfile); - - /* set up mutex and condition variable for i/o thread */ - /* hold mutex before creating i/o thread to keep it waiting. */ - pthread_mutex_init(&read_mutex, NULL); - pthread_mutex_init(&write_mutex, NULL); - pthread_cond_init(&write_cond, NULL); - - pthread_mutex_lock(&write_mutex); - buf_has_data=0; - pthread_mutex_unlock(&write_mutex); - - /* set up and launch i/o thread */ - pthread_attr_init(&iot_attr); - pthread_attr_setdetachstate(&iot_attr, PTHREAD_CREATE_JOINABLE); - pthread_create(&iothread, &iot_attr, &fix_imd_ioworker, this); - } -#endif -} - -/********************************* - * Clean up on deleting the fix. * - *********************************/ -FixIMD::~FixIMD() -{ - -#if defined(LAMMPS_ASYNC_IMD) - if (me == 0) { - pthread_mutex_lock(&write_mutex); - buf_has_data=-1; - pthread_cond_signal(&write_cond); - pthread_mutex_unlock(&write_mutex); - pthread_join(iothread, NULL); - - /* cleanup */ - pthread_attr_destroy(&iot_attr); - pthread_mutex_destroy(&write_mutex); - pthread_cond_destroy(&write_cond); - } -#endif - - inthash_t *hashtable = (inthash_t *)idmap; - memory->sfree(comm_buf); - memory->sfree(force_buf); - inthash_destroy(hashtable); - delete hashtable; - free(rev_idmap); - // close sockets - imdsock_shutdown(clientsock); - imdsock_destroy(clientsock); - imdsock_shutdown(localsock); - imdsock_destroy(localsock); - clientsock=NULL; - localsock=NULL; - return; -} - -/* ---------------------------------------------------------------------- */ -int FixIMD::setmask() -{ - int mask = 0; - mask |= POST_FORCE; - mask |= POST_FORCE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ -void FixIMD::init() -{ - if (strstr(update->integrate_style,"respa")) - nlevels_respa = ((Respa *) update->integrate)->nlevels; - - return; -} - -/* ---------------------------------------------------------------------- */ - -/* (re-)connect to an IMD client (e.g. VMD). return 1 if - new connection was made, 0 if not. */ -int FixIMD::reconnect() -{ - /* set up IMD communication, but only if needed. */ - imd_inactive = 0; - imd_terminate = 0; - - if (me == 0) { - if (clientsock) return 1; - if (screen && connect_msg) - if (nowait_flag) - fprintf(screen,"Listening for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); - else - fprintf(screen,"Waiting for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); - - connect_msg = 0; - clientsock = NULL; - if (nowait_flag) { - int retval = imdsock_selread(localsock,0); - if (retval > 0) { - clientsock = imdsock_accept(localsock); - } else { - imd_inactive = 1; - return 0; - } - } else { - int retval=0; - do { - retval = imdsock_selread(localsock, 60); - } while (retval <= 0); - clientsock = imdsock_accept(localsock); - } - - if (!imd_inactive && !clientsock) { - if (screen) - fprintf(screen, "IMD socket accept error. Dropping connection.\n"); - imd_terminate = 1; - return 0; - } else { - /* check endianness and IMD protocol version. */ - if (imd_handshake(clientsock)) { - if (screen) - fprintf(screen, "IMD handshake error. Dropping connection.\n"); - imdsock_destroy(clientsock); - imd_terminate = 1; - return 0; - } else { - int32 length; - if (imdsock_selread(clientsock, 1) != 1 || - imd_recv_header(clientsock, &length) != IMD_GO) { - if (screen) - fprintf(screen, "Incompatible IMD client version? Dropping connection.\n"); - imdsock_destroy(clientsock); - imd_terminate = 1; - return 0; - } else { - return 1; - } - } - } - } - return 0; -} - -/* ---------------------------------------------------------------------- */ -/* wait for IMD client (e.g. VMD) to respond, initialize communication - * buffers and collect tag/id maps. */ -void FixIMD::setup(int) -{ - /* nme: number of atoms in group on this MPI task - * nmax: max number of atoms in group across all MPI tasks - * nlocal: all local atoms - */ - int i,j; - int nmax,nme,nlocal; - int *mask = atom->mask; - int *tag = atom->tag; - nlocal = atom->nlocal; - nme=0; - for (i=0; i < nlocal; ++i) - if (mask[i] & groupbit) ++nme; - - MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); - maxbuf = nmax*size_one; - comm_buf = (void *) memory->smalloc(maxbuf,"imd:comm_buf"); - - connect_msg = 1; - reconnect(); - MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); - MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); - if (imd_terminate) - error->all("LAMMPS terminated on error in setting up IMD connection."); - - /* initialize and build hashtable. */ - inthash_t *hashtable=new inthash_t; - inthash_init(hashtable, num_coords); - idmap = (void *)hashtable; - - MPI_Status status; - MPI_Request request; - int tmp, ndata; - struct commdata *buf = static_cast(comm_buf); - - if (me == 0) { - int *taglist = new int[num_coords]; - int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ - - for (i=0; i < nlocal; ++i) { - if (mask[i] & groupbit) { - taglist[numtag] = tag[i]; - ++numtag; - } - } - - /* loop over procs to receive remote data */ - for (i=1; i < comm->nprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); - MPI_Send(&tmp, 0, MPI_INT, i, 0, world); - MPI_Wait(&request, &status); - MPI_Get_count(&status, MPI_BYTE, &ndata); - ndata /= size_one; - - for (j=0; j < ndata; ++j) { - taglist[numtag] = buf[j].tag; - ++numtag; - } - } - - /* sort list of tags by value to have consistently the - * same list when running in parallel and build hash table. */ - id_sort(taglist, 0, num_coords-1); - for (i=0; i < num_coords; ++i) { - inthash_insert(hashtable, taglist[i], i); - } - delete[] taglist; - - /* generate reverse index-to-tag map for communicating - * IMD forces back to the proper atoms */ - rev_idmap=inthash_keys(hashtable); - } else { - nme=0; - for (i=0; i < nlocal; ++i) { - if (mask[i] & groupbit) { - buf[nme].tag = tag[i]; - ++nme; - } - } - /* blocking receive to wait until it is our turn to send data. */ - MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status); - MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); - } - - return; -} - -/* worker threads for asynchronous i/o */ -#if defined(LAMMPS_ASYNC_IMD) -/* c bindings wrapper */ -void *fix_imd_ioworker(void *t) -{ - FixIMD *imd=(FixIMD *)t; - imd->ioworker(); - return NULL; -} - -/* the real i/o worker thread */ -void FixIMD::ioworker() -{ - while (1) { - pthread_mutex_lock(&write_mutex); - if (buf_has_data < 0) { - /* master told us to go away */ - fprintf(screen,"Asynchronous I/O thread is exiting.\n"); - buf_has_data=0; - pthread_mutex_unlock(&write_mutex); - pthread_exit(NULL); - } else if (buf_has_data > 0) { - /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); - } - delete[] msgdata; - buf_has_data=0; - pthread_mutex_unlock(&write_mutex); - } else { - /* nothing to write out yet. wait on condition. */ - pthread_cond_wait(&write_cond, &write_mutex); - pthread_mutex_unlock(&write_mutex); - } - } -} -#endif - -/* ---------------------------------------------------------------------- */ -/* Main IMD protocol handler: - * Send coodinates, energies, and add IMD forces to atoms. */ -void FixIMD::post_force(int vflag) -{ - /* check for reconnect */ - if (imd_inactive) { - reconnect(); - MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); - MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); - if (imd_terminate) - error->all("LAMMPS terminated on error in setting up IMD connection."); - if (imd_inactive) - return; /* IMD client has detached and not yet come back. do nothing. */ - } - - int *tag = atom->tag; - double **x = atom->x; - int *image = atom->image; - int nlocal = atom->nlocal; - int *mask = atom->mask; - struct commdata *buf; - - if (me == 0) { - /* process all pending incoming data. */ - int imd_paused=0; - while ((imdsock_selread(clientsock, 0) > 0) || imd_paused) { - /* if something requested to turn off IMD while paused get out */ - if (imd_inactive) break; - - int32 length; - int msg = imd_recv_header(clientsock, &length); - - switch(msg) { - - case IMD_GO: - if (screen) - fprintf(screen, "Ignoring unexpected IMD_GO message.\n"); - break; - - case IMD_IOERROR: - if (screen) - fprintf(screen, "IMD connection lost.\n"); - /* fallthrough */ - - case IMD_DISCONNECT: { - /* disconnect from client. wait for new connection. */ - imd_paused = 0; - imd_forces = 0; - memory->destroy(force_buf); - force_buf = NULL; - imdsock_destroy(clientsock); - clientsock = NULL; - if (screen) - fprintf(screen, "IMD client detached. LAMMPS run continues.\n"); - - connect_msg = 1; - reconnect(); - if (imd_terminate) imd_inactive = 1; - break; - } - - case IMD_KILL: - /* stop the simulation job and shutdown IMD */ - if (screen) - fprintf(screen, "IMD client requested termination of run.\n"); - imd_inactive = 1; - imd_terminate = 1; - imd_paused = 0; - imdsock_destroy(clientsock); - clientsock = NULL; - break; - - case IMD_PAUSE: - /* pause the running simulation. wait for second IMD_PAUSE to continue. */ - if (imd_paused) { - if (screen) - fprintf(screen, "Continuing run on IMD client request.\n"); - imd_paused = 0; - } else { - if (screen) - fprintf(screen, "Pausing run on IMD client request.\n"); - imd_paused = 1; - } - break; - - case IMD_TRATE: - /* change the IMD transmission data rate */ - if (length > 0) - imd_trate = length; - if (screen) - fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); - break; - - case IMD_ENERGIES: { - IMDEnergies dummy_energies; - imd_recv_energies(clientsock, &dummy_energies); - break; - } - - case IMD_FCOORDS: { - float *dummy_coords = new float[3*length]; - imd_recv_fcoords(clientsock, length, dummy_coords); - delete[] dummy_coords; - break; - } - - case IMD_MDCOMM: { - int32 *imd_tags = new int32[length]; - float *imd_fdat = new float[3*length]; - imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); - - if (imd_forces < length) { /* grow holding space for forces, if needed. */ - memory->destroy(force_buf); - force_buf = (void *) memory->smalloc(length*size_one, - "imd:force_buf"); - } - imd_forces = length; - buf = static_cast(force_buf); - - /* compare data to hash table */ - for (int ii=0; ii < length; ++ii) { - buf[ii].tag = rev_idmap[imd_tags[ii]]; - buf[ii].x = imd_fdat[3*ii]; - buf[ii].y = imd_fdat[3*ii+1]; - buf[ii].z = imd_fdat[3*ii+2]; - } - delete[] imd_tags; - delete[] imd_fdat; - break; - } - - default: - if (screen) - fprintf(screen, "Unhandled incoming IMD message #%d. length=%d\n", msg, length); - break; - } - } - } - - /* update all tasks with current settings. */ - int old_imd_forces = imd_forces; - MPI_Bcast(&imd_trate, 1, MPI_INT, 0, world); - MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); - MPI_Bcast(&imd_forces, 1, MPI_INT, 0, world); - MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); - if (imd_terminate) - error->all("LAMMPS terminated on IMD request."); - - if (imd_forces > 0) { - /* check if we need to readjust the forces comm buffer on the receiving nodes. */ - if (me != 0) { - if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ - if (force_buf != NULL) - memory->sfree(force_buf); - force_buf = memory->smalloc(imd_forces*size_one, "imd:force_buf"); - } - } - MPI_Bcast(force_buf, imd_forces*size_one, MPI_BYTE, 0, world); - } - - /* Check if we need to communicate coordinates to the client. - * Tuning imd_trate allows to keep the overhead for IMD low - * at the expense of a more jumpy display. Rather than using - * end_of_step() we do everything here in one go. - * - * If we don't communicate, only check if we have forces - * stored away and apply them. */ - if (update->ntimestep % imd_trate) { - if (imd_forces > 0) { - double **f = atom->f; - buf = static_cast(force_buf); - - /* XXX. this is in principle O(N**2) == not good. - * however we assume for now that the number of atoms - * that we manipulate via IMD will be small compared - * to the total system size, so we don't hurt too much. */ - for (int j=0; j < imd_forces; ++j) { - for (int i=0; i < nlocal; ++i) { - if (mask[i] & groupbit) { - if (buf[j].tag == tag[i]) { - f[i][0] += imd_fscale*buf[j].x; - f[i][1] += imd_fscale*buf[j].y; - f[i][2] += imd_fscale*buf[j].z; - } - } - } - } - } - return; - } - - /* check and potentially grow local communication buffers. */ - int i, k, nmax, nme=0; - for (i=0; i < nlocal; ++i) - if (mask[i] & groupbit) ++nme; - - MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); - if (nmax*size_one > maxbuf) { - memory->sfree(comm_buf); - maxbuf = nmax*size_one; - comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); - } - - MPI_Status status; - MPI_Request request; - int tmp, ndata; - buf = static_cast(comm_buf); - - if (me == 0) { - /* collect data into new array. we bypass the IMD API to save - * us one extra copy of the data. */ - msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; - msgdata = new char[msglen]; - imd_fill_header((IMDheader *)msgdata, IMD_FCOORDS, num_coords); - /* array pointer, to the offset where we receive the coordinates. */ - float *recvcoord = (float *) (msgdata+IMDHEADERSIZE); - - /* add local data */ - if (unwrap_flag) { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double xy = domain->xy; - double xz = domain->xz; - double yz = domain->yz; - - for (i=0; i> 10 & 1023) - 512; - int iz = (image[i] >> 20) - 512; - - if (domain->triclinic) { - recvcoord[j] = x[i][0] + ix * xprd + iy * xy + iz * xz; - recvcoord[j+1] = x[i][1] + iy * yprd + iz * yz; - recvcoord[j+2] = x[i][2] + iz * zprd; - } else { - recvcoord[j] = x[i][0] + ix * xprd; - recvcoord[j+1] = x[i][1] + iy * yprd; - recvcoord[j+2] = x[i][2] + iz * zprd; - } - } - } - } - } else { - for (i=0; inprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); - MPI_Send(&tmp, 0, MPI_INT, i, 0, world); - MPI_Wait(&request, &status); - MPI_Get_count(&status, MPI_BYTE, &ndata); - ndata /= size_one; - - for (k=0; kxprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double xy = domain->xy; - double xz = domain->xz; - double yz = domain->yz; - - for (i=0; i> 10 & 1023) - 512; - int iz = (image[i] >> 20) - 512; - - if (domain->triclinic) { - buf[nme].tag = tag[i]; - buf[nme].x = x[i][0] + ix * xprd + iy * xy + iz * xz; - buf[nme].y = x[i][1] + iy * yprd + iz * yz; - buf[nme].z = x[i][2] + iz * zprd; - } else { - buf[nme].tag = tag[i]; - buf[nme].x = x[i][0] + ix * xprd; - buf[nme].y = x[i][1] + iy * yprd; - buf[nme].z = x[i][2] + iz * zprd; - } - ++nme; - } - } - } else { - for (i=0; i(num_coords+maxbuf+imd_forces)*size_one; -} - -/* End of FixIMD class implementation. */ - -/***************************************************************************/ - -/* NOTE: the following code is the based on the example implementation - * of the IMD protocol API from VMD and NAMD. The UIUC license allows - * to re-use up to 10% of a project's code to be used in other software */ - -/*************************************************************************** - * DESCRIPTION: - * Socket interface, abstracts machine dependent APIs/routines. - ***************************************************************************/ - -int imdsock_init(void) { -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) - int rc = 0; - static int initialized=0; - - if (!initialized) { - WSADATA wsdata; - rc = WSAStartup(MAKEWORD(1,1), &wsdata); - if (rc == 0) - initialized = 1; - } - - return rc; -#else - return 0; -#endif -} - - -void * imdsock_create(void) { - imdsocket * s; - - s = (imdsocket *) malloc(sizeof(imdsocket)); - if (s != NULL) - memset(s, 0, sizeof(imdsocket)); - - if ((s->sd = socket(PF_INET, SOCK_STREAM, 0)) == -1) { - printf("Failed to open socket."); - free(s); - return NULL; - } - - return (void *) s; -} - -int imdsock_bind(void * v, int port) { - imdsocket *s = (imdsocket *) v; - memset(&(s->addr), 0, sizeof(s->addr)); - s->addr.sin_family = PF_INET; - s->addr.sin_port = htons(port); - - return bind(s->sd, (struct sockaddr *) &s->addr, sizeof(s->addr)); -} - -int imdsock_listen(void * v) { - imdsocket *s = (imdsocket *) v; - return listen(s->sd, 5); -} - -void *imdsock_accept(void * v) { - int rc; - imdsocket *new_s = NULL, *s = (imdsocket *) v; -#if defined(ARCH_AIX5) || defined(ARCH_AIX5_64) || defined(ARCH_AIX6_64) - unsigned int len; -#define _SOCKLEN_TYPE unsigned int -#elif defined(SOCKLEN_T) - SOCKLEN_T len; -#define _SOCKLEN_TYPE SOCKLEN_T -#elif defined(_POSIX_SOURCE) - socklen_t len; -#define _SOCKLEN_TYPE socklen_t -#else -#define _SOCKLEN_TYPE int - int len; -#endif - - len = sizeof(s->addr); - rc = accept(s->sd, (struct sockaddr *) &s->addr, ( _SOCKLEN_TYPE * ) &len); - if (rc >= 0) { - new_s = (imdsocket *) malloc(sizeof(imdsocket)); - if (new_s != NULL) { - *new_s = *s; - new_s->sd = rc; - } - } - return (void *)new_s; -} - -int imdsock_write(void * v, const void *buf, int len) { - imdsocket *s = (imdsocket *) v; -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) - return send(s->sd, (const char*) buf, len, 0); /* windows lacks the write() call */ -#else - return write(s->sd, buf, len); -#endif -} - -int imdsock_read(void * v, void *buf, int len) { - imdsocket *s = (imdsocket *) v; -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) - return recv(s->sd, (char*) buf, len, 0); /* windows lacks the read() call */ -#else - return read(s->sd, buf, len); -#endif - -} - -void imdsock_shutdown(void *v) { - imdsocket * s = (imdsocket *) v; - if (s == NULL) - return; - -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) - shutdown(s->sd, SD_SEND); -#else - shutdown(s->sd, 1); /* complete sends and send FIN */ -#endif -} - -void imdsock_destroy(void * v) { - imdsocket * s = (imdsocket *) v; - if (s == NULL) - return; - -#if defined(_MSC_VER) || defined(__MINGW32_VERSION) - closesocket(s->sd); -#else - close(s->sd); -#endif - free(s); -} - -int imdsock_selread(void *v, int sec) { - imdsocket *s = (imdsocket *)v; - fd_set rfd; - struct timeval tv; - int rc; - - if (v == NULL) return 0; - - FD_ZERO(&rfd); - FD_SET(s->sd, &rfd); - memset((void *)&tv, 0, sizeof(struct timeval)); - tv.tv_sec = sec; - do { - rc = select(s->sd+1, &rfd, NULL, NULL, &tv); - } while (rc < 0 && errno == EINTR); - return rc; - -} - -int imdsock_selwrite(void *v, int sec) { - imdsocket *s = (imdsocket *)v; - fd_set wfd; - struct timeval tv; - int rc; - - if (v == NULL) return 0; - - FD_ZERO(&wfd); - FD_SET(s->sd, &wfd); - memset((void *)&tv, 0, sizeof(struct timeval)); - tv.tv_sec = sec; - do { - rc = select(s->sd + 1, NULL, &wfd, NULL, &tv); - } while (rc < 0 && errno == EINTR); - return rc; -} - -/* end of socket code. */ -/*************************************************************************/ - -/*************************************************************************/ -/* start of imd API code. */ -/* Only works with aligned 4-byte quantities, will cause a bus error */ -/* on some platforms if used on unaligned data. */ -void swap4_aligned(void *v, long ndata) { - int *data = (int *) v; - long i; - int *N; - for (i=0; i>24)&0xff) | ((*N&0xff)<<24) | - ((*N>>8)&0xff00) | ((*N&0xff00)<<8)); - } -} - - -/** structure used to perform byte swapping operations */ -typedef union { - int32 i; - struct { - unsigned int highest : 8; - unsigned int high : 8; - unsigned int low : 8; - unsigned int lowest : 8; - } b; -} netint; - -static int32 imd_htonl(int32 h) { - netint n; - n.b.highest = h >> 24; - n.b.high = h >> 16; - n.b.low = h >> 8; - n.b.lowest = h; - return n.i; -} - -static int32 imd_ntohl(int32 n) { - netint u; - u.i = n; - return (u.b.highest << 24 | u.b.high << 16 | u.b.low << 8 | u.b.lowest); -} - -static void imd_fill_header(IMDheader *header, IMDType type, int32 length) { - header->type = imd_htonl((int32)type); - header->length = imd_htonl(length); -} - -static void swap_header(IMDheader *header) { - header->type = imd_ntohl(header->type); - header->length= imd_ntohl(header->length); -} - -static int32 imd_readn(void *s, char *ptr, int32 n) { - int32 nleft; - int32 nread; - - nleft = n; - while (nleft > 0) { - if ((nread = imdsock_read(s, ptr, nleft)) < 0) { - if (errno == EINTR) - nread = 0; /* and call read() again */ - else - return -1; - } else if (nread == 0) - break; /* EOF */ - nleft -= nread; - ptr += nread; - } - return n-nleft; -} - -static int32 imd_writen(void *s, const char *ptr, int32 n) { - int32 nleft; - int32 nwritten; - - nleft = n; - while (nleft > 0) { - if ((nwritten = imdsock_write(s, ptr, nleft)) <= 0) { - if (errno == EINTR) - nwritten = 0; - else - return -1; - } - nleft -= nwritten; - ptr += nwritten; - } - return n; -} - -int imd_disconnect(void *s) { - IMDheader header; - imd_fill_header(&header, IMD_DISCONNECT, 0); - return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); -} - -int imd_handshake(void *s) { - IMDheader header; - imd_fill_header(&header, IMD_HANDSHAKE, 1); - header.length = IMDVERSION; /* Not byteswapped! */ - return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); -} - -/* The IMD receive functions */ - -IMDType imd_recv_header(void *s, int32 *length) { - IMDheader header; - if (imd_readn(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE) - return IMD_IOERROR; - swap_header(&header); - *length = header.length; - return IMDType(header.type); -} - -int imd_recv_mdcomm(void *s, int32 n, int32 *indices, float *forces) { - if (imd_readn(s, (char *)indices, 4*n) != 4*n) return 1; - if (imd_readn(s, (char *)forces, 12*n) != 12*n) return 1; - return 0; -} - -int imd_recv_energies(void *s, IMDEnergies *energies) { - return (imd_readn(s, (char *)energies, sizeof(IMDEnergies)) - != sizeof(IMDEnergies)); -} - -int imd_recv_fcoords(void *s, int32 n, float *coords) { - return (imd_readn(s, (char *)coords, 12*n) != 12*n); -} - -// Local Variables: -// mode: c++ -// compile-command: "make -j4 openmpi" -// c-basic-offset: 2 -// fill-column: 76 -// indent-tabs-mode: nil -// End: - diff --git a/src/USER-IMD/fix_imd.h b/src/USER-IMD/fix_imd.h deleted file mode 100644 index a7c3cb36eb..0000000000 --- a/src/USER-IMD/fix_imd.h +++ /dev/null @@ -1,131 +0,0 @@ -/* -*- 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - The FixIMD class contains code from VMD and NAMD which is copyrighted - by the Board of Trustees of the University of Illinois and is free to - use with LAMMPS according to point 2 of the UIUC license (10% clause): - -" Licensee may, at its own expense, create and freely distribute -complimentary works that interoperate with the Software, directing others to -the TCBG server to license and obtain the Software itself. Licensee may, at -its own expense, modify the Software to make derivative works. Except as -explicitly provided below, this License shall apply to any derivative work -as it does to the original Software distributed by Illinois. Any derivative -work should be clearly marked and renamed to notify users that it is a -modified version and not the original Software distributed by Illinois. -Licensee agrees to reproduce the copyright notice and other proprietary -markings on any derivative work and to include in the documentation of such -work the acknowledgement: - - "This software includes code developed by the Theoretical and Computational - Biophysics Group in the Beckman Institute for Advanced Science and - Technology at the University of Illinois at Urbana-Champaign." - -Licensee may redistribute without restriction works with up to 1/2 of their -non-comment source code derived from at most 1/10 of the non-comment source -code developed by Illinois and contained in the Software, provided that the -above directions for notice and acknowledgement are observed. Any other -distribution of the Software or any derivative work requires a separate -license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to -negotiate an appropriate license for such distribution." -------------------------------------------------------------------------- */ - -#ifdef FIX_CLASS - -FixStyle(imd,FixIMD) - -#else - -#ifndef LMP_FIX_IMD_H -#define LMP_FIX_IMD_H - -#include "fix.h" - -#if defined(LAMMPS_ASYNC_IMD) -#include -#endif - -/* prototype for c wrapper that calls the real worker */ -extern "C" void *fix_imd_ioworker(void *); - -namespace LAMMPS_NS { - -class FixIMD : public Fix { - public: - FixIMD(class LAMMPS *, int, char **); - ~FixIMD(); - int setmask(); - void init(); - void setup(int); - void post_force(int); - void post_force_respa(int, int, int); - double memory_usage(); - - protected: - int imd_port; - void *localsock; - void *clientsock; - - int num_coords; // total number of atoms controlled by this fix - int size_one; // bytes per atom in communication buffer. - int maxbuf; // size of atom communication buffer. - void *comm_buf; // communication buffer - void *idmap; // hash for mapping atom indices to consistent order. - int *rev_idmap; // list of the hash keys for reverse mapping. - - int imd_forces; // number of forces communicated via IMD. - void *force_buf; // force data buffer - double imd_fscale; // scale factor for forces. in case VMD's units are off. - - int imd_inactive; // true if IMD connection stopped. - int imd_terminate; // true if IMD requests termination of run. - int imd_trate; // IMD transmission rate. - - int unwrap_flag; // true if coordinates need to be unwrapped before sending - int nowait_flag; // true if LAMMPS should not wait with the execution for VMD. - int connect_msg; // flag to indicate whether a "listen for connection message" is needed. - - int me; // my MPI rank in this "world". - int nlevels_respa; // flag to determine respa levels. - - int msglen; - char *msgdata; - -#if defined(LAMMPS_ASYNC_IMD) - int buf_has_data; // flag to indicate to the i/o thread what to do. - pthread_mutex_t write_mutex; // mutex for sending coordinates to i/o thread - pthread_cond_t write_cond; // conditional variable for coordinate i/o - pthread_mutex_t read_mutex; // mutex for accessing data receieved by i/o thread - pthread_t iothread; // thread id for i/o thread. - pthread_attr_t iot_attr; // i/o thread attributes. -public: - void ioworker(void); -#endif - -protected: - int reconnect(); -}; - -} - -#endif -#endif - -// Local Variables: -// mode: c++ -// compile-command: "make -j4 openmpi" -// c-basic-offset: 2 -// fill-column: 76 -// indent-tabs-mode: nil -// End: diff --git a/src/USER-MISC/README b/src/USER-MISC/README index e868ed4ac4..d3724d17b4 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -4,7 +4,7 @@ file (actually a *.cpp and *.h file). More information about each feature can be found by reading its doc page in the LAMMPS doc directory. This link points to the doc -page for all LAMMPS commands: +page for all LAMMPS input script commands: http://lammps.sandia.gov/doc/Section_commands.html#3_5\ @@ -20,8 +20,12 @@ angle_style cosine/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 angle_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 +compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 +fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009 +fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008 pair dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11 +pair eam/cd, Alexander Stukowski, stukowski at mm.tu-darmstadt.de, pair lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 diff --git a/src/USER-ACKLAND/compute_ackland_atom.cpp b/src/USER-MISC/compute_ackland_atom.cpp similarity index 100% rename from src/USER-ACKLAND/compute_ackland_atom.cpp rename to src/USER-MISC/compute_ackland_atom.cpp diff --git a/src/USER-ACKLAND/compute_ackland_atom.h b/src/USER-MISC/compute_ackland_atom.h similarity index 100% rename from src/USER-ACKLAND/compute_ackland_atom.h rename to src/USER-MISC/compute_ackland_atom.h diff --git a/src/USER-SMD/Install.sh b/src/USER-SMD/Install.sh deleted file mode 100644 index 954dda9102..0000000000 --- a/src/USER-SMD/Install.sh +++ /dev/null @@ -1,15 +0,0 @@ -# Install/unInstall package files in LAMMPS - -if (test $1 = 1) then - - cp -p fix_smd.cpp .. - - cp -p fix_smd.h .. - -elif (test $1 = 0) then - - rm -f ../fix_smd.cpp - - rm -f ../fix_smd.h - -fi diff --git a/src/USER-SMD/README b/src/USER-SMD/README deleted file mode 100644 index d4dbc8513d..0000000000 --- a/src/USER-SMD/README +++ /dev/null @@ -1,20 +0,0 @@ -The files in this directory are a user-contributed package for LAMMPS. - -The person who created these files is Axel Kohlmeyer -(akohlmey@cmm.chem.upenn.edu). Contact him directly if you have -questions or for example scripts that use it. - -This package implements a "fix smd" command which can be used in a -LAMMPS input script. SMD stands for steered molecular dynamics. - -Here is what Axel said about the new fix: - -Please find attached a fix that I was asked to write for our -"coarse-grainers". They need free energy profiles for their -parametrization and this seems to work best. The code basically is a -version of "fix spring" on steroids (constant velocity pulling with -velocity set to zero should recover "fix spring"), but I found it -cleaner to have it in a seprarate class since otherwise the result may -become too convoluted for people that only want to use "fix spring". - - diff --git a/src/USER-SMD/fix_smd.cpp b/src/USER-SMD/fix_smd.cpp deleted file mode 100644 index 4d251d3fe3..0000000000 --- a/src/USER-SMD/fix_smd.cpp +++ /dev/null @@ -1,403 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Axel Kohlmeyer (UPenn) - based on fix spring by: Paul Crozier (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdlib.h" -#include "string.h" -#include "fix_smd.h" -#include "atom.h" -#include "comm.h" -#include "update.h" -#include "respa.h" -#include "domain.h" -#include "error.h" -#include "group.h" - -using namespace LAMMPS_NS; - -enum { SMD_NONE=0, - SMD_TETHER=1<<0, SMD_COUPLE=1<<1, - SMD_CVEL=1<<2, SMD_CFOR=1<<3, - SMD_AUTOX=1<<4, SMD_AUTOY=1<<5, SMD_AUTOZ=1<<6}; - -#define SMALL 0.001 - -/* ---------------------------------------------------------------------- */ - -FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - styleflag = SMD_NONE; - k_smd = f_smd = v_smd = -1.0; - xflag = yflag = zflag = 1; - xc = yc = zc = 0.0; - xn = yn = zn = 1.0; - pmf = r_old = r_now = r0 = 0.0; - - restart_global = 1; - vector_flag = 1; - size_vector = 7; - global_freq = 1; - extvector = 1; - - int argoffs=3; - if (strcmp(arg[argoffs],"cvel") == 0) { - if (narg < argoffs+3) error->all("Illegal fix smd command"); - styleflag |= SMD_CVEL; - k_smd = atof(arg[argoffs+1]); - v_smd = atof(arg[argoffs+2]); // to be multiplied by update->dt when used. - argoffs += 3; - } else if (strcmp(arg[argoffs],"cfor") == 0) { - if (narg < argoffs+2) error->all("Illegal fix smd command"); - styleflag |= SMD_CFOR; - f_smd = atof(arg[argoffs+1]); - argoffs += 2; - } else error->all("Illegal fix smd command"); - - if (strcmp(arg[argoffs],"tether") == 0) { - if (narg < argoffs+5) error->all("Illegal fix smd command"); - styleflag |= SMD_TETHER; - if (strcmp(arg[argoffs+1],"NULL") == 0) xflag = 0; - else xc = atof(arg[argoffs+1]); - if (strcmp(arg[argoffs+2],"NULL") == 0) yflag = 0; - else yc = atof(arg[argoffs+2]); - if (strcmp(arg[argoffs+3],"NULL") == 0) zflag = 0; - else zc = atof(arg[argoffs+3]); - r0 = atof(arg[argoffs+4]); - if (r0 < 0) error->all("R0 < 0 for fix smd command"); - argoffs += 5; - } else if (strcmp(arg[argoffs],"couple") == 0) { - if (narg < argoffs+6) error->all("Illegal fix smd command"); - styleflag |= SMD_COUPLE; - igroup2 = group->find(arg[argoffs+1]); - if (igroup2 == -1) - error->all("Could not find fix smd couple group ID"); - if (igroup2 == igroup) - error->all("Two groups cannot be the same in fix smd couple"); - group2bit = group->bitmask[igroup2]; - - if (strcmp(arg[argoffs+2],"NULL") == 0) xflag = 0; - else if (strcmp(arg[argoffs+2],"auto") == 0) styleflag |= SMD_AUTOX; - else xc = atof(arg[argoffs+2]); - if (strcmp(arg[argoffs+3],"NULL") == 0) yflag = 0; - else if (strcmp(arg[argoffs+3],"auto") == 0) styleflag |= SMD_AUTOY; - else yc = atof(arg[argoffs+3]); - if (strcmp(arg[argoffs+4],"NULL") == 0) zflag = 0; - else if (strcmp(arg[argoffs+4],"auto") == 0) styleflag |= SMD_AUTOZ; - else zc = atof(arg[argoffs+4]); - - r0 = atof(arg[argoffs+5]); - if (r0 < 0) error->all("R0 < 0 for fix smd command"); - argoffs +=6; - } else error->all("Illegal fix smd command"); - - force_flag = 0; - ftotal[0] = ftotal[1] = ftotal[2] = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int FixSMD::setmask() -{ - int mask = 0; - mask |= POST_FORCE; - mask |= POST_FORCE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::init() -{ - double xcm[3], xcm2[3]; - masstotal = group->mass(igroup); - group->xcm(igroup,masstotal,xcm); - - double dx,dy,dz; - if (styleflag & SMD_TETHER) { - dx = xc - xcm[0]; - dy = yc - xcm[1]; - dz = zc - xcm[2]; - } else { /* SMD_COUPLE */ - masstotal2 = group->mass(igroup2); - group->xcm(igroup2,masstotal2,xcm2); - if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; - else dx = xc; - if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; - else dy = yc; - if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; - else dz = zc; - } - - if (!xflag) dx = 0.0; - if (!yflag) dy = 0.0; - if (!zflag) dz = 0.0; - r_old = sqrt(dx*dx + dy*dy + dz*dz); - if (r_old > SMALL) { - xn = dx/r_old; - yn = dy/r_old; - zn = dz/r_old; - } - - if (strstr(update->integrate_style,"respa")) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::setup(int vflag) -{ - if (strstr(update->integrate_style,"verlet")) - post_force(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::post_force(int vflag) -{ - if (styleflag & SMD_TETHER) smd_tether(); - else smd_couple(); - - if (styleflag & SMD_CVEL) r_old += v_smd * update->dt; -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::smd_tether() -{ - double xcm[3]; - group->xcm(igroup,masstotal,xcm); - - // fx,fy,fz = components of k * (r-r0) - - double dx,dy,dz,fx,fy,fz,r,dr; - - dx = xcm[0] - xc; - dy = xcm[1] - yc; - dz = xcm[2] - zc; - r_now = sqrt(dx*dx + dy*dy + dz*dz); - - if (!xflag) dx = 0.0; - if (!yflag) dy = 0.0; - if (!zflag) dz = 0.0; - r = sqrt(dx*dx + dy*dy + dz*dz); - if (styleflag & SMD_CVEL) { - if(r > SMALL) { - double fsign; - fsign = (v_smd<0.0) ? -1.0 : 1.0; - dr = r - r0 - r_old; - fx = k_smd*dx*dr/r; - fy = k_smd*dy*dr/r; - fz = k_smd*dz*dr/r; - pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt; - } - } else { - r_old = r; - fx = f_smd*dx/r; - fy = f_smd*dy/r; - fz = f_smd*dz/r; - } - - // apply restoring force to atoms in group - // f = -k*(r-r0)*mass/masstotal - - double **f = atom->f; - int *mask = atom->mask; - int *type = atom->type; - double *mass = atom->mass; - int nlocal = atom->nlocal; - - ftotal[0] = ftotal[1] = ftotal[2] = 0.0; - force_flag = 0; - - double massfrac; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massfrac = mass[type[i]]/masstotal; - f[i][0] -= fx*massfrac; - f[i][1] -= fy*massfrac; - f[i][2] -= fz*massfrac; - ftotal[0] -= fx*massfrac; - ftotal[1] -= fy*massfrac; - ftotal[2] -= fz*massfrac; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::smd_couple() -{ - double xcm[3],xcm2[3]; - group->xcm(igroup,masstotal,xcm); - group->xcm(igroup2,masstotal2,xcm2); - - // renormalize direction of spring - double dx,dy,dz,r,dr; - if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; - else dx = xn*r_old; - if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; - else dy = yn*r_old; - if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; - else dz = zn*r_old; - if (!xflag) dx = 0.0; - if (!yflag) dy = 0.0; - if (!zflag) dz = 0.0; - r = sqrt(dx*dx + dy*dy + dz*dz); - if (r > SMALL) { - xn = dx/r; yn = dy/r; zn = dz/r; - } - - double fx,fy,fz; - if (styleflag & SMD_CVEL) { - dx = xcm2[0] - xcm[0]; - dy = xcm2[1] - xcm[1]; - dz = xcm2[2] - xcm[2]; - r_now = sqrt(dx*dx + dy*dy + dz*dz); - - dx -= xn*r_old; - dy -= yn*r_old; - dz -= zn*r_old; - - if (!xflag) dx = 0.0; - if (!yflag) dy = 0.0; - if (!zflag) dz = 0.0; - r = sqrt(dx*dx + dy*dy + dz*dz); - dr = r - r0; - - if (r > SMALL) { - double fsign; - fsign = (v_smd<0.0) ? -1.0 : 1.0; - - fx = k_smd*dx*dr/r; - fy = k_smd*dy*dr/r; - fz = k_smd*dz*dr/r; - pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt; - } - - } else { - dx = xcm2[0] - xcm[0]; - dy = xcm2[1] - xcm[1]; - dz = xcm2[2] - xcm[2]; - r_now = sqrt(dx*dx + dy*dy + dz*dz); - r_old = r; - - fx = f_smd*xn; - fy = f_smd*yn; - fz = f_smd*zn; - } - - // apply restoring force to atoms in group - // f = -k*(r-r0)*mass/masstotal - - double **f = atom->f; - int *mask = atom->mask; - int *type = atom->type; - double *mass = atom->mass; - int nlocal = atom->nlocal; - - ftotal[0] = ftotal[1] = ftotal[2] = 0.0; - force_flag = 0; - - double massfrac; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - massfrac = mass[type[i]]/masstotal; - f[i][0] += fx*massfrac; - f[i][1] += fy*massfrac; - f[i][2] += fz*massfrac; - ftotal[0] += fx*massfrac; - ftotal[1] += fy*massfrac; - ftotal[2] += fz*massfrac; - } - if (mask[i] & group2bit) { - massfrac = mass[type[i]]/masstotal2; - f[i][0] -= fx*massfrac; - f[i][1] -= fy*massfrac; - f[i][2] -= fz*massfrac; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::write_restart(FILE *fp) -{ -#define RESTART_ITEMS 5 - double buf[RESTART_ITEMS], fsign; - - if (comm->me == 0) { - // make sure we project the force into the direction of the pulling. - fsign = (v_smd<0.0) ? -1.0 : 1.0; - buf[0] = r_old; - buf[1] = xn*fsign; - buf[2] = yn*fsign; - buf[3] = zn*fsign; - buf[4] = pmf; - int size = RESTART_ITEMS*sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(&buf[0],sizeof(double),RESTART_ITEMS,fp); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::restart(char *buf) -{ - double *list = (double *)buf; - r_old = list[0]; - xn=list[1]; - yn=list[2]; - zn=list[3]; - pmf=list[4]; -} - -/* ---------------------------------------------------------------------- */ - -void FixSMD::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- - return components of total smd force on fix group -------------------------------------------------------------------------- */ - -double FixSMD::compute_vector(int n) -{ - // only sum across procs one time - - if (force_flag == 0) { - MPI_Allreduce(ftotal,ftotal_all,3,MPI_DOUBLE,MPI_SUM,world); - force_flag = 1; - if (styleflag & SMD_CVEL) { - ftotal_all[3]=ftotal_all[0]*xn+ftotal_all[1]*yn+ftotal_all[2]*zn; - ftotal_all[4]=r_old; - } else { - ftotal_all[3]=f_smd; - ftotal_all[4]=r_old; - } - ftotal_all[5]=r_now; - ftotal_all[6]=pmf; - } - return ftotal_all[n]; -} diff --git a/src/USER-SMD/fix_smd.h b/src/USER-SMD/fix_smd.h deleted file mode 100644 index 5319ca9558..0000000000 --- a/src/USER-SMD/fix_smd.h +++ /dev/null @@ -1,60 +0,0 @@ -/* ---------------------------------------------------------------------- - 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(smd,FixSMD) - -#else - -#ifndef LMP_FIX_SMD_H -#define LMP_FIX_SMD_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixSMD : public Fix { - public: - FixSMD(class LAMMPS *, int, char **); - int setmask(); - void init(); - void setup(int); - void post_force(int); - void post_force_respa(int, int, int); - double compute_vector(int); - - void write_restart(FILE *); - void restart(char *); - - private: - double xc,yc,zc,xn,yn,zn,r0; - double k_smd,f_smd,v_smd; - int xflag,yflag,zflag; - int styleflag; - double r_old,r_now,pmf; - - int igroup2,group2bit; - double masstotal,masstotal2; - int nlevels_respa; - double ftotal[3],ftotal_all[7]; - int force_flag; - - void smd_tether(); - void smd_couple(); -}; - -} - -#endif -#endif